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Preface 


This report comprises the final report for a study- 
entitled, "Transform Processing and Coding of Images, n performed 
by the Electronic Sciences Laboratory of the University of Southern 
California for the Jet Propulsion Laboratory under JPL Contract 
952312. Mr. Thomas Rindfleisch of JPL served as project director 
for the study. This report supplants the interim report USCEE No. 
341 entitled '‘Transform Processing and Coding of Images, * 
published in March, 1969. Pertinent introductory material from 
that report is included in the present report for completeness. 
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1. Introduction 


The basic goal of digital image coding is the development 
of a coding technique that permits the representation, and 
subsequent recovery, of an image by a minimal number of code 
bits [l-3]. In some applications virtually no image distortion is 
permitted in the coding process, while in other applications a 
controlled amount of distortion is allowable in the achievement of 
a substantial bit reduction. In general, when redundancy is 
removed from a data source, the compressed data is more 
sensitive to the effect of channel errors. One of the restrictions 
in selecting a coding method, therefore, is that the compressed 
data must not be overly sensitive to channel errors. 

In 1967 a new technique of image coding, called Fourier 
transform^ coding, was developed at the University of Southern 
California [4-6], Another related method, called Hadamard transform 
coding, was discovered at USC in 1968 [7-8], Since then investigations 
have been made into the applications of other mathematical 
transforms for image coding. Out of these studies has emerged 
the generalized technique of transform image coding [-9-11]. 
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1. 1 Image Transform Coding 


Figure 1-1 contains a block diagram of the image 
transform coding system. In operation a two-dimensional 
transform is taken of the brightness samples of an image, or 
subsection of an image, on a line by line basis. The resultant 
transform samples are then operated upon by a sample selector 
that selects which samples are to be transmitted on the basis 
of magnitude or position in the plane. Those samples that are 
to be transmitted are quantized and coded. At the receiver, the 
data is decoded to reconstruct the transform domain, and an 
inverse transform is taken to reconstruct the original image. 

A bandwidth reduction is achieved simply by not trans- 
mitting all of the transform domain samples. Those samples 
that are not transmitted are generally of such low magnitude that 
they contribute little in the image reconstruction. 

There are two basic forms of sample selection- -zonal 
sampling and thre shold. sampling --that can be employed. In 
zonal sampling, only those transform samples that lie within 
certain geometric region in- the 'transform domain are selected 
for transmission. The basic problem with zonal sampling is that 
in certain pictures many large magnitude samples may lie without 
the zonal region and will, therefore, not be transmitted. In order 
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Figure 1-1. Generalized Transform Coding of Images 










to avoid such errors it is possible to establish a threshold level 
on the magnitude of transform domain samples such that if the 
transform sample magnitude is greater than the threshold it 
will be selected, and the sample will be deleted if it falls below 
the threshold. With threshold coding it is necessary to code the 
location in the transform domain of a selected sample as well 
as its value. 

The major advantage of image transform coding other than 
its potential for bandwidth compression is the tolerance to channel 
errors that transform coding affords. An intuitive justification 
for transmitting the transform of an image rather than the spatial 
representation of the image is that for many transforms the 
channel noise introduced in the image transform tends to be 
distributed evenly over the entire reconstructed image. Consequently, 
the channel noise is manifested as a low spatial frequency error 
in reconstruction. Experimental evidence indicates that the eye 
is more sensitive -to the high frequency discrete errors caused 
by channel errors in the spatial domain than it is to the same 
number of errors in the transform domain. 
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1,2 Original Images 


Figure 1-2 contains photographs o£ the five original 
images that have been used as test images for the evaluation 
of image transform coding. These images contain 256 by 256 
elements quantized to 64 grey levels. The images were read 
from magnetic tape, displayed on a Hewlett-Packard Model 1300 
cathode ray tube display, and photographed with Polaroid Type 
47 film. 


- 4 ' 



a * Surveyor footpad 


b. Moonscape 




c. Surveyor experimental box 


d« Surveyor boom 



e , Girl 


Figure 1-2 Original Test Images 
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Z« Image Transformation 


In this section consideration is given to the mathematical 
formulation of image transforms. The characteristics and 
properties of the Fourier, Hadamard, Haar, Karhunen-L-oeve, and 
a class of transitional transforms are briefly developed. Experi- 
mental results are presented. 

2. 1 Formulation 

An image may be represented by an array of intensity 

components or samples over the image surface by two dimensional 

sampling. For the present discussion an image array will be 

Z 

considered to be a square array of N intensity samples described 
by the function f(x, y) over the image coordinates (Xpy). 

Conceptually, there are two major types of image transforms 
which shall be called transforms of the first and second kind, A 
transform of the first kind maps a two dimensional image array 

of dimension N X N into a one dimensional vector of dimension 

2 

1 X N according to the relation 
N-l N-l 

F(w) = 2 2 f(x, y) a(x, y, w) <2~I) 

X— U y— 0 
2 

for w = 0, 1, 2, * * - , N -I 
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where a[x,y,w) is the forward transform kernel of the first 
kind. A reverse transform of the first kind is defined as 


Nf-3 

l(x,y) = £i F(wJ b(x, y, w) (Z- Z) 

w= 0 

for x, y = 0, 1, 2, - * * , N— 1 

where b(x p y,w] is the reverse transform kernel of the first 
kind. A transform of the second kind maps an image array of 
dimension N X N into a two dimensional array of the same 
dimension as defined by 

N - 1 N- 1 

F(u,v) = 2 lJ H f(x, y] a(x, y, u# v) (2-3) 

x~0 y=0 

for u, v = 0, 1 , Z t “ * ■ * N — 1 

where a(x, y t u f v) is the forward transform of the second kind. 

The corresponding reverse transformation is given by 

A N-l N- 1 

f(x # y) - 1/ li F(u,v) b(x, y p u,v) (2-4) 

u= 0 v— 0 

for x, y - 0, 1, N— I 

where b(x t y,u p v) is the reverse transform kernel of the second 
kind. For transforms of the first and second kind, when the 
function f (x p y) resulting from the reverse transform operation is 
equivalent to the original image, f(x p y) t the reverse transform is 
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called an inverse transform. Transforms of the first and second 
kind are said to be orthogonal if the following conditions 
are met* 


Transforms of the first kind: 


S a(x,y t w) a* (Qf, B, w) = 6(x-a p y-B) 
w 

(2- 5a) 

b(x, y,w) b*(£* p B,w) = y-B) 

w 

(2- 5b) 

S £a(x, y, w) a* <x P y, 9) = 6 (w-G) 

X y 

<2-5c) 

^EMx. YpW) tf*(x p y, 9) = 6(w-^ 

Xy 

(2- 5d) 

Transforms of the second kind; 


2 2 afx.yjUpv) a*(a p B, u p v) = 6(x-af,y-B) 
u v 

*Z-6a) 

2 2 b{x p y j, Up v) b+fa, B, u p v) - 6 (x-a p y-B) 
u v 

(2- 6b) 

2 2 a(x p y $ Up v) a*(x, y t 9 P Cp) = &(u-8p v-q>) 

X y 

(2- 6c) 

2 2 b(x p y„ Up v) h*{xp y r @ P q>) = &(u-9»v-cp) 

x y 

{2- 6d) 


The limits of summation are eliminated in subsequent 
equations unless required for clarity. 
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A forward transform kernel of the second kind is said 


to be separable if it can be written as 

a(x,y,u f v) = a^Y, v) (21-7) 

A separable two dimensional transform can be computed in two 
steps. First, a one dimensional transform is taken along each 
row of the image, f(x,y) p yielding 

N- 1 

F(u,y)=2 f(x,y)a(x,u) (2-8) 

x=0 1 

Next, a second one dimensional transform is taken along each 
column of F(u, y) giving 

N - 1 

F(u,v) = £ F(u, y) a,(y, v) (2-9) 

y=0 L 

The transformation kernel is called separable symmetric if 

y,u,v) = a^x.u) a (y f v) (2-10) 

For ease of implementation, the separable symmetric property is 
desirable. 

It is often useful to express two dimensional transforms in 
matrix notation. For example, with a forward transform kernel of 
the second kind that is separable symmetric let: 
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[£] - image matrix, f(x, y) 

[f] = transformed image matrix, F(u, v) 

[a] = transform matrix, A(o , P ) 

Then by matrix multiplication 

CF] = [a] [£] [A] (2-11) 

Now pre- and post- multiplication of each side of [F] by a 
reverse transform matrix, [b], gives 

fi] ^ [B] [F] [B] - [B] [A][f][A] [B] (2-12) 

where [I] is, in general, an approximation of &]♦ If the reverse 
transform matrix is the inverse matrix [A] ] of [A], then 

[f] = [A]' 1 [AKf] CaKA]' 1 (2-13) 

But 

[A]' 1 [A] = [A] [A]* 1 = [I] (2-14) 

where £l] is the identify matrix* Hence 

[£]= It} = [A]' 1 CfHa]' 1 (2- IB) 

Thus, f(x, y) and F(u, v) can be expressed as two dimensional 
transform pairs if [a] has an inverse. If [A] is a unitary matrix, 
then by definition 


- 10 - 


unitary matrix 


(2-16} 


[B] 5 [A]’ 1 = [ A ]* t 

where [A]* is the complex conjugate matrix of [a] and [a] T 
is the matrix transpose of [A]. If in addition [A] is symmetric 

[b] = [a] = [A] symmetric unitary matrix (2-17) 

A real, unitary matrix is called an orthogonal matrix* For 
such a matrix 


m - [ a ] -1 

Finally, if [a] 
[B> [A]' 1 


T 

= LA] orthogonal matrix 

is a symmetric orthogonal matrix, then 

- [A] symmetric orthogonal matrix 


( 2 - 18 ) 


(2-19) 


If the forward transformation matrix is constrained to be 
unitary, then the transformation can be interpreted as a decomposition 
of the image data into a generalized two dimensional spectrum. 

Each spectral component in the transform domain corresponds to 
the amount of energy of the spectral function within the original 
image. In this context the concept of frequency may now be 
generalized to include transformations of functions other than sine 
and cosine waveforms. This type of generalized spectral analysis 
is useful in the investigation of specific decompositions which are 
best suited for particular classes of images. 
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The following paragraphs contain an analysis of the Fourier 


Hadamard, Haar, transitional, and Karhunen-Loeve transformations 
with particular emphasis on their applicability to image processing. 

2, 2 Fourier Transform 

The discrete Fourier transform with and without efficient 
computational algorithms, has long been used for signal 
analysis [12]. Only recently have Fourier transform methods 
been utilized for image coding [4-6]. 

The two dimensional Fourier transform of an image field, 
f(x, y), maybe expressed as 


The inverse Fourier transform which reconstructs the original 
image is given by 


Since the transform kernels are separable and symmetric the two 
dimensional transform can be computed as two sequential one 
dimensional transforms. 


i N-l N-l 

F(u» v) = — S S f(x, y) exp 

N x =0 y-0 



) 


( 2 - 20 ) 



( 2 - 21 ) 
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The terms u and v are called the spatial frequencies 
of the image in analogy with time series analysis. When the 
Fourier transform relationship is expressed in the form given 
by Equation (2-20) the origin, or zero spatial frequency term appears 
in the corner of the transform plane. For display purposes it is 
convenient to shift the origin to the center of the transform 
domain. This is easily accomplished by multiplying the image by 
the function (-l) X+y before the transformation [13]. 

Even though f(x, y) is a real positive function, its transform, 

F(u,v), is in general complex. Thus, while the image contains 

2 2 
N components, the transform contains 2N components* the real 

and imaginary, or magnitude and phase components of each 

spatial frequency. However, since f{x, y) is a real positive function, 

F(u,v) exhibits a property of conjugate symmetry [l3]. Specifically, 

F(u,v) = F*( -u, -v) (2-22) 

As a result of the conjugate symmetry property of the Fourier trans- 
form it is only necessary to transmit the samples of one half of the 
transform plane; the other half can be reconstructed from the half plane 

3ft 

samples transmitted * Hence* the Fourier transform of an image 


A reconstruction of the original can be obtained from the half 

plane transform samples directly by a Hilbert filtering technique [13], 
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can be described by N data components. 

The two dimensional Fourier transform of an image is 
essentially a Fourier series representation of a two dimensional 
field. For the Fourier series representation to be valid the field 
must be periodic. Thus, the original image must be considered to 
be periodic horizontally and vertically. The right side of the image 
therefore abuts the left side and the top and bottom of the image are 
adjacent. Spatial frequencies along the coordinate axes of the 
transform plane arise from, these transitions. Although these are 
false spatial frequencies from the standpoint of being necessary for 
representing the image within the image boundary, they do not impair 
reconstruction. On the contrary, these spatial frequencies are 
required to reconstruct the sharp boundaries of the image. 

Figure 2-1 presents displays of the Fourier transforms 
in shifted form of two of the original test scenes. The logarithm 
of the magnitude of each transform is displayed rather than the 
magnitude itself in order to reduce the dynamic range of the display. 
In addition, a threshold display is presented in which all the absolute 
values above the threshold are set to white and all others are made 
black. Such a display gives a graphic illustration of the heavy 
concentration of energy around the origin (center of photograph) of 
the Fourier transform. 
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a. Logarithm of the magnitude of the 
Surveyor box transform 


b* Threshold display of the Surveyor 
box transform 




c* Logarithm of the magnitude of 
tne moonscape transform 



d. Threshold display of the moon- 
scape transform 


Figure 2-1 Fourier Transforms 
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2, 3 Hadamard Transform 


The Hadamard transform, also known as the Walsh 
transform, is based upon the Hadamard matrix which is a square 
array of plus and minus ones whose rows and columns are orthogonal 
to one another [ 14-16] - If [ H] is an N by N Hadamard matrix, 
then the product of N and its transpose is 

Ch][H] T =N[I] (2-23) 

If [H]is a symmetric Hadamard matrix, then Equation (2-23) 
reduces to 

[H] [H] = N[l] (2-24) 


A Hadamard matrix multiplied by the normalization factor 


1 


is an orthonormal matrix. 


The lowest order Hadamard matrix is the Hadamard matrix 



(2-25) 


It is known that if a Hadamard matrix of order N exists (N >2), then 
N = 0 (mod 4). The existence of a Hadamard matrix for every 
value of N satisfying this requirement has not been shown, but 
constructions are available for nearly all permissible values of N 
up to 200, The simplest construction is for a Hadamard matrix 




of order N s £ n where n ia an integer* In this case if [H^l 
is a Hadamard matrix of order N, the matrix 


H 


N 


H 


N 



H, 


-H. 


(Z-26) 


N 


N — 1 


ia a Hadamard matrix of order 2N. 

A frequency interpretation can be given to the Hadamard 

matrix generated from the core matrix of Equation (2-25)* Along 

each row of the Hadamard matrix the frequency is called the number 

of changes in sign* Harmuth has coined the word ,r sequency M to 

designate the number of sign changes [17] . It is possible to 

construct a Hadamard matrix of order N = Z n that has frequency 

components at every integer from 0 to N-I* 

This frequency interpretation of the rows of a Hadamard 

matrix leads one to consider the rows to be equivalent to 

rectangular waves ranging between ±1 with a sub-period of — 

N 

units* Such functions are called Walsh functions [18-22] and are 
fur die r related to the Rademacher functions [23] * Thus, in 
this context the Hadamard matrix merely performs the decomposition 
of a function by a set of rectangular waveforms rather than the 
sine-cosine waveforms associated with the Fourier transform* 

For symmetric Hadamard matrices of order N = 2°, the 
two dimensional Hadamard transform may be written in series 
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form as 


T N-l N-I 

Fiu,v)= ±2 E f(x, y)(-l) px ’ y ’ u ' v) 
^ x=0 y=0 


( 2 - 27 ) 


n-l 

where p(x, y, u, v) = . (u.x. + v.y.), The term a o., v.,x.. and y. 

are the binary representationa of u, v, ,x, and y respectively. 

For example. 


1U) DECIMAL = ( Vl “n-2 *" U 1 “oWnARY 


where e {0 # l] . 

Another aeries representation exists for a Hadamard matrix 
in "ordered" form in which the sequency of each row is larger 
than the preceding row. By this representation 


F(u,v) = 


ir 

N x=0 


E ’f(x.y) M> q(X - V ' U,V) 

y=0 


(Z-Z8) 


where 

n-l 

q(x, y, u ( v) = £ [g.(u) x, + g.(v) yj 

i=0 21 1 

and 


g 0 <u> 

= 

u 1 
n-l 

gjfu) 

= 

n - 1 n- 2 

g 2 {«} 

s 

u - + u * 
n-2 n-3 

W u) 

m 

u + u 

1 0 
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The two dimensional Hadamard transform may be computed 
in either natural or ordered form with an algorithm analogous 
to the fast Fourier transform computer algorithm* 

Figure 2-2 presents the ordered Hadamard transforms 
of two test scenes. The origin of the transform domain is now in 
the lower left corner and the axes are now spatial sequencies as 
opposed to spatial frequencies. Notice that as in the Fourier 
case, the image energy tends to concentrate itself heavily in the 
lowest spatial sequency areas providing the potential for large 
band w idth r e due tion s , Aga in both a lo g a r ithmi c an d thr e s ho Id 
display are provided for dynamic range purposes, 

2, 4 Transitional Transforms 

In related work [24-27] it has been shown that a class of 
rapidly implementable orthogonal transformations exists for matrices 
composed of Kronecker products of smaller core matrices. In 
fact, both the Hadamard and Fourier transforms have been shown 
to be subsets of this much larger class of Kronecker transform 
matrices* A class of transformations exists for which the Hadamard 
transform is a limiting case, and these transforms will now be 
investigated as to their image processing potential. 

The transformation resulting from performing the Kronecker 
operation of the core matrix 
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a . Logarithm of the magnitude of the 
Surveyor box transform 


b. Threshold display of the Surveyor 
box transform 










c. Logarithm of the magnitude of 
the moonscape transform 


d. Threshold display of the moon- 
scape transform 


Figure 2-2 Hadamard Transforms 
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[H] = 


cos 8 sin8 
sin 3 -cos 8 


] 


(2-29} 


with itself n times results in a matrix whose row and column 
entries, indexed by x and u respectively, can be described 
by the following equation 


of the column and row indexes respectively* It is evident that 
while the Hadamard transform has received considerable attention 
(often under the name of the discrete Walsh transform) it is 
important to note that this transform is the limiting case of the 
powers of two Krone eke r transforms presented above* As 8 

O Q 

varies between 0 and 45 the transforms vary from a diagonal 

matrix to the Hadamard matrix at 45°. In the process of varying 0 

over this interval, the transformations have ranged from having 

o 

all of their energy on the diagonal at 0 to uniform energy spread 
at 45° (Hadamard case). Figure 2-3 presents examples of the 
transitional transforms of the Surveyor box test scene for four 
different values of 0, Notice that functions of the magnitude are 
displayed in all cases because even for the 8 = 0° diagonal case, 
negative signs on the diagonal are possible. It is evident from 



(2-30) 


The u^ and variables are the bits in the binary representation 
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H01 



c« 0 - 30°, threshold display 
max. value = 6,611 


cos © 
sin © 



d. © = 45°, Hadamard transform, 
threshold display, max. 
value = 11,486 

sin 0 
-cos © 



Figure 2-3 Transitional Transforms 
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this figure that the transform which computes the image into 
the fewest significant coefficients is the Hadamard transform* 


Z. S Haar Transform 


The Haar transform [28] ia another transformation that] 
like the Hadamard or Walsh transform, requires no multiplications* 
The Haar matrix consists of plus and minus ones as well as zeros 
and is non- symmetric, orthogonal but not orthonormal (unless 
multiplied by the proper diagonal matrix). The Haar matrix can be 
likened to a sampling system in which various rows sample the 
input with finer and finer resolution increasing in powers of two* 

An 8x8 orthonormal Haar matrix is shown below: 


l 

1 

1 

1 

l 

l 

1 

1 

1 

1 

1 

1 

-I 

-1 

-1 

-1 

41 

Jz 

~4z 

-42 

0 

0 

0 

0 

0 

0 

0 

0 

41 

42 

V 2 

-41 

2 

-2 

0 

0 

0 

0 

0 

0 

0 

0 

2 

-2 

0 

0 

0 

0 

0 

0 

0 

0 

2 

-2 

0 

0 

0 

0 

0 

0 

0 

0 

2 

-2 


The Haar transform is defined for data of resolution equal 
to a power of two, and the matrix is factorable into a product of 
matrices with a large number of zero entries [25], Consequently, 
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a fast algorithm also exists for this transform. The number of 
computer operations required for a vector matrix multiplication 
is given by as compared to the 2N log^N requirement 

of Fourier or Hadamard, which itself is a considerable savings 
over the normal vector -matrix multiplication requirement of N Z 
operations. As with the Walsh functions, the Haar functions 
can be generalized to contain entries of roots of unity other 
than ±1, W atari [Z9] has described the generalized Haar system 
and has shown that it is possible to preserve some of the original 
Haar convergence properties. The extension to matrix factorization 
is straightforward and will not be pursued further. However, the 
number of operations necessary to implement a p th order generalized 
Haar transform is given by a geometric progression resulting in 
p(N-l) /(p-1), In image processing applications, the Haar transform 
provides a transform domain in which a type of differential energy 
is concentrated in localized regions. Thus there is an area in 
which adjacent picture element differential energy is concentrated, 
(the upper right quarter of the transform plane), an area in which 
differential energy of adjacent picture elements taken two at a 
time is concentrated, and in general an area in which difference 
energy of adjacent picture elements taken a power of two at a time 
is concentrated. 
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Figure 2-4 presents the Haar transform of the test 
scenes. The logarithmic results vividly display the derivative 
energy effect especially in the upper right quarter of the plane. 
Note that in the Haar transform there is also a concentration of 
image energy in the lower left corner or origin of the transform 
plane. The data point at the origin in the Fourier p Hadamard, and 
Haar transforms all are equal to the average energy in the original 
image and correspond to the row of all r, ones ,T in the transform 
matrices, 

2,6 Karhunen-Loeve Transform 

The Karhunen-Loeve transform is a special case of an 
eigenvector matrix transformation [30-37], Consider a real 
symmetric matrix [c] of order n. The eigenvectors of [c] 
are column vectors [K.] p i = 1* 2, - ■ - f n satisfying the relationship 

[ c3 [k.3 = x. [kJ 

i ii 

where the scalars X . are the eigenvalues of [c] , Let a square 
matrix [K] p called the modal matrix of [c] p he constructed 
from the eigenvector columns in the following manner: 

[K] = [tKj] [K z ] ■ • • tK n ]] 


Also let the eigenvalues be located along the diagonal of a matrix 



a * Logarithm of the magnitude of the 
Surveyor box transform 


b. Threshold display of the Surveyor 
box transform 




c. Logarithm of the magnitude of 
the moonscape transform 


d. Threshold display of the moon- 
scape transform 


Figure 2-4 Haar Transforms 
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Then by equation (2-32) 

[C][K] = [ K] [e] 

Now, pre multiplication of equation (2-34) by 1 k 3 1 gives 
[K] _1 [C][K] = [K]' 1 [K][E] = [E] 

Taking the transpose of both sides of equation (2-35) yields 

[k] t [c] t [tK] _1 ] T = [e] t 

But, since [C] is a symmetric matrix and [E] is diagonal, by 
correspondence between equation (2-35) and (2-36) 

W 1 = [K] T 

Thus, if they exist, the eigenvectors of a matrix are orthogonal. 
It can be easily shown [38] that when [c] is symmetric, its 
eigenvalues are all real quantities. 

Consider now a data column vector [f] of length m. The 
eigenvector transform [f] of [f] is then 


(2-34) 


(2-35) 


(2-36) 


(2-37) 
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[F] = [K] [f] 


(2-38) 


and the inverse eigenvector transform [f ] of [f] is 

[£ ]= [K] T [F] = [K] T = [£] (2-39) 

Thus [f] and [F] are transform pairs of an orthogonal matrix 
transformation. The vector [F] represents a matrix 
decomposition of [f ] into a set of orthogonal waveforms defined 
by [K] . Generally, the exact form of the orthogonal functions 
cannot be easily described. 

If only the first q of the m columns of [K] are employed in 
the forward and reverse transform, then the mean square error 
between the original and the reconstructed data vector is [30, 31 ] 

m 

S = 2 (2-40) 

k=q+l k 

Since the are monotonically decreasing in value, the error 
will be minimum for any q. 

When the eigenvector matrix [K] is composed of 
eigenvectors of the covariance matrix 

[ C] « E | [f(i)-«ir ] C f(j)-7(jj] | (2-41) 

for i, j = 1, 2, • • • , n, of the data vectors, the resulting 

eigenvector matrix [K] is called the Karhunen- Loeve (K-L) transform. 
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For a two dimensional image transform of the first kind* the 
forward transform kernel a(x* y, w) of Equation (2-1) satisfies 
the equation 


M-l N-l 

X£w) a(x,y,w) = £ Tt C[x,x‘ * y* y 1 ] a(x' f y , ,wl 
x'= 0 y'-0 


(2-42) 



ponding K-L transform of the second kind described by the kernel 
a(x* y, Uj v) of Equation (2-3) is found from 


for u, v = 0, 1 B - ** , N— 1, where X (u* v) are a two dimensional 
ordering of the eigenvalues X(w), If the covariance function in 
Equation (2-43) can be written as 


then the transform kernel a(x, y, u ( v) cam be separated. The 
resulting two dimensional transform can then be computed sequentially 
along each row and column of the image. 

Figure 2-5 contains photographs of an image that has been 
Karhunen-Loeve transformed in 4X4 element blocks. The 16x16 


N-1N-1 

X(ujv) a(x p y*u,v) = £ £ ctx.x 1 , y,y'} a(x’ t y% u, t) 

x’=0 y'=0 


■ (2-43) 


Cfxjx'jy.y^s C 1 fx p x'} C^ty.y 1 } 


(2-44) 
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a. Logarithm of the 
magnitude of the 
Surveyor box transform 


b. Inverse transform of 
transform of Surveyor box 




c. Logarithm of the 
magnitude of the 
girl transform 

A - 1 . 

B = .8 

G ” . 6 

D - E = F = 0 


d* Inverse transform 
of transform of girl 


Figure 2-5 Karhunen-Loeve Transforms in 4 x A Element Blocks 
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component correlation matrix from which fee K-L transform 
matrix was derived is shown in Figure Z-6 . 

2 . 7 Summary 

For image coding the desirable properties of a mathematical 
transform are that the transform redistribute the image energy 
to as few transform domain samples as possible, and furthermore 
feat the transform be easily computable. The Fourier and 
Hadamard transforms fulfill both requirements, and will be 
analyzed in greater detail in subsequence sections. 

None of the transitional transforms, other than fee 
Hadamard transform, provide a compact distribution of energy. 
These transforms will not be considered further for image coding, 
but it should be noted feat the transitional transforms may be 
useful for dimensionality reduction for pattern recognition 
applications. 

The Haar transform possesses an extremely fast 
computational algorithm. However, the peculiar spatial sampling 
procedure-- sampling in pairs of elements --does not appear 
to be particularly useful for image coding, and therefore fee Haar 
transform will not be considered further. The Haar transform 
may find some usefulness, however, for digital edge 
enhancement since fee transform domain is a mapping of the spatial 
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12 3 4 


5 6 7 8 

9 10 11 12 

13 14 15 16 

a * Element Array 


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

1 {a B C D B C D E' C 5 I F D E F G~J 

2 | B ABCC BGDDCDEEDEF 

3jCBA B D CB CE DC DF E DE 

4 S D C B A E DCBFEDCGF ED 

s;bcdea BCDBCDECDEF 
6CBCDB ABCC BCDDC DE 

7 DCB CC BABDCBCE DCD 

8 E DCB D CBA E DCBF E DC 

9 C D E F B CDEABCDBCDE 

10 - D C D E C B CDBA B CCB CD 

11 ' E D ODD CB CCBA BDC BC 

12FEDCE DCBDCBAEDCB 

13 DE FGC DE FBCDEA B CD 

14 E DE FDCDE CBCDBA BC 

15FEDFE DCDDGB CCB AB 

16 G F F D F FDCEDCBDCBAi 

l 1 

b. Correlation Matrix 


Figure 2-6 Correlation Matrix Model for 4 x 4 Element 
Karhunen-Loeve Transform 
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differential energy of the original image. 

The Karhunen-Loeve transform provides the best 
compaction of image energy for natural images. The major 
difficulty associated with the use of the Karhunen-Loeve 
transform for image coding is the great amount of computation 
involved. First, the image correlation matrix must be 

estimated or modeled. Next, the correlation matrix must be diagonalized 
to determine its eigenvalues and eigenvectors. Finally, the transform 
itself must be taken. In general » there is no fast computational 
algorithm for the transform. In those applications in which the 
amount of computation is not of principal concern, the Karhunen-Loeve 
transform may find practical application. Furthermore, since the 
K-L transform is the optimum image transform in a mean square 
error sense, when sample deletion is employed, it is worthwhile 
to consider its performance as a standard for other image 
transforms. 

The next two sections contain a general analysis of the 
Fourier, Hadamard, and Karhunen-Loeve image transforms. 
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3. Statistical Analysis of Image Transforms 

The development of efficient quantization and coding methods 
for image transform samples requires an understanding of the 
statistical properties of the transform domain. This section 
presents a derivation of the first and second moments of transform 
samples, and also contains the development of a stochastic model 
for the probability density of transform samples. 

The statistical analysis of image transforms is predicated 
on the representation of an original image as a two dimensional 
stochastic process, f(x, y). The spatial mean 

E{f(x,y)} ■ f(x,y) {3-1) 

and the covariance 

E |[f(x 1 ,y 1 )-f(x 1 ,y 1 )] [f(x 2 ,y 2 )-f(x 2 ,y 2 ) ■ Cfa^.x^ y^ y 2 J 

(3-2) 

are assumed known or at least estimateable. Appendix A 
describes measurements of the covariance function of an 
image. 

3. 1 Moments 

Kor a generalized forward transform given by 
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(3-3) 


N-l N-l 

F(u,v) = ^ f(x,y) a(x, y,u,v) 

the mean of the transform samples is simply the forward 
transform of the mean of the image samples. Thus, 

N-l N-l 

E{F(u,v)}= F(u,v) - Jj Jj f(x, y) a(x, y, u, v) 
x=0 y=0 

For an ordered, orthonormal transform with an average value 
term 


F(u,v) = Nf(x,y) 6 (u, v) 


The covariance function of the transform domain samples 
is by definition 


C{ V u 2’ V 1 >V 2 } = E 



U l )-F( V V l ,][F * (ll 2’ 


v 2 )- f <u 2 , 



Substitution of equations (3-3) and (3-4) yields 


C{u 1 ,u 2 ,v 1 ,v 2 } = E S t£(x 1 ,y 1 )-f(x 1 ,y 1 )] a(x 1> y 1 ,u 1 ,v 1 )J 

&*(x 2> y 2 )-f(x 2> y 2 )*] a^(x 2 ,y 2 ,u 2 ,v 2 )]| 


or 


(3-4a) 


(3- 4b) 


(3-5) 


(3-6) 


- 35 - 



(3-7) 


c{u ,U ,v ,v }= D E E E E [f(x ,y )-f(x ,y )] 

1 1 2 X 1 71*2*2 l 1 1 11 

• [r*(x 2J y 2 )-f(x 2 ,y 2 r]j a(x 1 > y 1 ,u 1 ,v 1 ) a*(x 2> y^, v 2 ) 

The expected value of the bracketed term'in the. summation of 
Equation (3-7) is by definition the spatial domain covariance function, 
Cfxj, x^, y , y 2 ). Hence, 

C{u , u 2 > v i» v o) = s SS S C{x ,x ,y ,y ) 

x i x 2 y i y 2 

•■(Xj.yj.ttj.Vj) a*(x 2 ,y 2 ,u 2> v 2 ) (3-8) 


The variance of the transform domain samples is 
2 

a (u, v) = c[u,u,v, v} (3-9) 

Therefore, the general expression for the variance of transform 
domain samples becomes 


a 2 (u,v)= SEED C{x ,x ,y ,y } 

x i x 2 yi^ 

•a(x 1 ,y 1 ,u, v) a*(x 2 ,y 2 ,u, v) (3-10) 

There are two special cases, of interest. For an image that is 
statistically stationary in the spatial domain 
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(3-11) 


Ctx l’ x 2’ y i’ y 2 3 = C t x 1 ~ X 2’ y i" y 2 3 

If the original image is uncorrelated in the horizontal and 
vertical directions 

C[x 1 ,x 2> y 1 ,y 2 } = C^fxj, x 2 } c (3-12) 

2 

and the transform domain variance can be computed as 0 (u, v) = 

2 2 

CT (u) o (v) provided that the transform kernel is separable. 

Further investigation of the variance of transform domain samples 
requires specification of the transform. 

Fourier Transform 

For the Fourier transform the variance function of Equation 
(3-10) can be written as 


£ '{u J v)= £ 2 2 £ C{x lf x 2 , y lt y 2 } exp / - [u(x 1 -x 2 ) 

l y 2 V 

+ v( yi -y 2 )]l 


N Xj x 2 Tl y 2 


(3-13) 


Consider the case for which the original image is stationary 
and orthogonally uncorrelated. The variance function reduces to 


2 2 2 
o (u, v) = a (u) 0 (v) 


(3-14) 


where 
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(3-15) 


a 2 (p) = -Jf p S C{ z r z 2 } exp / - ^ p ( Z r » 2 )\ 

Z 1 Z 2 * 

with p = XL or v and z. = x. or y... The coordinate variance 
111 

function may then be rewritten as 


N 


C 2 (P) = jj S exp ~ P z 9 ^ £ C{ Zl -z,}exp^ - ^ p Zl 


1 2 


N 


The second summation is the one dimensional discrete Fourier 
transform of the covariance function shifted by z^. By the 
Fourier transform translation theorem 


° 2 (P) 




2TTi 

N 



G(p) 


or 


0 2 (p) = G(p) 


where G(p) and C{ ! z^} are one dimensional discrete Fourier 

transform pairs. If the transform is over a complete image 

dimension, G(p) is the discrete version of the power spectral 

density, ^Mp), of the image function along one coordinate minus 

the average image power, S (0), along the coordinate. Hence, 

z 

C 2 (P) = S (p)— S (0) 

Z Z 


(3-16) 


(3- 17a) 


(3- 17b) 


(3-18) 
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Thus, the transform domain sample variance along a coordinate 
direction is directly proportional to the power spectral density 
of the image along the corresponding orthogonal coordinate. 

If the original image function can be considered to be a Gauss - 
Markov process, the covariance function is [39] 

C z (z l' z 2 ) = C z <°> exp | _Y I z r z 2 l | (3-19) 

where C z (0) is a scaling constant and y is a shape constant. 

Then 

S z (p)-S z (0) = C (0) r -fi— 1 (3-20) 

L Y + p J 

For the Markov process example, the transform domain variance 
becomes 


2 , 

a (u, v) 


<yo) c y (0) 


4 g 3 

(a 2 + u 2 )(P 2 +v 2 ) 


(3-21) 


where *-^(0) a-^d C y {0) are the magnitude scaling constants 

and the shape constants of the spatial domain covariance function, 

respectively. 


Hadamard transform 

From Equations (3-10) and (2-20) the variance function 
of transform domain samples for the ordered Hadamard transform is 
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a 2 (u,v)= -i- £ £ £ £ Cfxj.x^y .y J 

x i x 2 y i y 2 

“- 1 r -1 

i= S 0 L g i (U)(X li + X 2i } + g i (v)(y ii + (3- 22) 

(- 1 ) 


Since the Hadamard transform does not possess a sequency shifting 
property it is not possible to reduce Equation (3-22) to closed' 
form. 

Karh.unen-L.oeve Transform 

The general expression for the variance of the variance 
of transform samples given by Equation (3-10) can be rewritten as 


a (u,v) = £ D a*(x ,y u,v)£ £ C{x ,x y y } a(x , y , u, v) 

x 2 y 2 *1^! 


(3-23) 


For the Karhunen-Loeve transform from Equation (2-43) the second 
set of summations defines the transform kernel. Thus, 


\(9,cp) a(x 2 , y 2 , 0,cp) = 


D E C{x 1 ,x 2 ,y 1 ,y 2 } a(x ,y ,0,(p) 
x i 


(3-24) 


where X(0.,cp) are the eigenvalues of the covariance matrix. By 
this equivalence 
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(3-25) 



Since the Karhunen-Loeve transform is an orthogonal 
transformation, from Equation (2- 6c) 

a 2 (u, v) = X (u, v) (3-26) 

and the variance of each transform sample is equal to its 
corresponding eigenvalue. 

3. 2 Probability Densities 

It would be desirable to know the probability density . of 
transform samples for an arbitrary image transform. Unfortunately, 
this result is not easily obtained since the original image probability 
density is not usually well defined, and also, the transform 
operation is quite often mathematically complex. However, the 
transforms considered for image processing applications form a 
weighted sum over all of the elements in the original image. 

Therefore, one can evoke qualitative arguments based upon the 
Central Limit. Theorem of statistics that the probability density of 
transform samples tends to be Gaussian with moments as calculated 
in the previous section,. For the subsequent analysis, a Gaussian 
model is developed for the probability density of the Fourier, 

Hadamard, and Karhunen-Loeve transform domain samples. 
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Fourier transform samples are complex numbers which 
may be represented in real and imaginary, or magnitude and phase, 
form. In. either case there are two components per transform 
sample that must be quantized. The real, F R (u,v) and imaginary, 
Fj(u, v), components of the Fourier transform samples may be 

assumed to follow the same Gaussian distribution whose variance, 

2 ’ 

.O' (u, v), is proportional to the power spectral density of the original 

image. Hence, 


p \ F r (u,v) 


p \ Fj(u, v) 


O 1 f (U,v) ^ 

[2Tia (u,v)] ^ exp 


[ 2tt0^(u-, v)] ^ exp ^ 



(3-27) 


(3-28) 


If the real and imaginary components are Gaussian, the magnitude 
of the Fourier transform sample, F^Xu, v), is Rayleigh distributed 


P 



2 , 

o (u, v) 


exp 


- f m K v) 

2o 2 (u, v)- 


V u * v) > 0 


(3-29) 


audits phase, F R (u, v), is uniformly distributed 
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p 




2tt 


— tt < Fp < -f tt 


( 3 - 30 ) 


Hadamard transform samples are real, bipolar numbers which can 

be represented by a single component per sample. The statistical 

distribution of Hadamard sample components, F„(u, v), maybe 

H 

considered to follow a Gaussian distribution of the form 


P 


'H 



[ 2tt a (u, v) ] 2 exp 


-F h (u,v) 
I 2a 2 (u,v) 


( 3 - 31 ) 


Karhunen-Loeve transform samples are also real bipolar numbers. 
The probability density of the samples maybe modeled as 


P 


K 



[ 2tt O (u, v) ] ^ exp 


/ _f k ^ 

|2a 2 (u, v) 


( 3 - 32 ) 


When the variance function, C (u, v), is not known for a 
particular image, or class of images,, to be transformed, the 
function can usually .be modeled without seriously affecting the 
quantization process. From examination of the Fourier, Hadamard, 
and Karhunen-Hoeve transforms of a typical image, it can be 
deduced that the variance function should be a maximum at the 
origin in the transform domain, be circularly symmetric, and 
decrease in magnitude monotonically toward the higher spatial 
frequencies. A two dimensional function processing these 
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characteristics is the Gaussian shaped curve described by 


2 2 
O (u, v) = S exp ^ 


2 2 
u -f v 

p/2 


where S is an amplitude scaling constant and p is a spread 
control constant. Another useful function- for modeling of the 
variance function is 


0 (u, v) = 


2 2 2 2 
(u + a ) (v + 0 ) 


where S is an amplitude scaling constant and a and 0 are 
spread control constants. This model holds exactly for the 
Fourier transform if the original image can be considered as a 
Gauss -Markov process source. 


( 3 - 33 ) 


( 3 - 34 ) 
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4. Generalized Transform Coding 


The basic premise of image transform coding is that the 
two dimensional transform of an image has an energy distribution 
more amenable to coding than the spatial domain representation. 

As a result of the inherent element- to- element correlation of 
natural images, for many image transforms, the energy in the 
transform domain tends to be clustered in a relatively few number 
of transform samples. This property can be exploited to achieve 
a sample reduction compared to conventional spatial domain coding. 

There are two methods of obtaining a sample reduction by 
transform coding- -zonal sampling and threshold sampling. In 
zonal sampling the image reconstruction is made with a subset, 
usually the lowest spatial coefficients, of the- transform domain 
samples. Those samples which are employed in the reconstruction 
are chosen before the transformation on the basis of expected energy. 
With threshold sampling the reconstruction is made with a subset 
of the largest magnitude transform domain samples. 

This section presents a discussion of the performance of 
the Karhunen-Loeve, Fourier, and Hadamard transforms for zonal 
and threshold sampling in the transform domain. The three transforms 
are compared on the basis of minimum mean square error. Experi- 
mental results are presented for a subjective comparison. 
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4. 1 Generalized Zonal Sampling 
Optimum Zonal Sampling 


Consider an image transform of the first kind. With 

zonal sampling the image is reconstructed with the first M of 
2 

the N transform samples. Thus, the reconstructed image becomes 
M- 1 

f(x, y) = E F(w)b(x,y,w) (4-1) 

w=0 

The mean square error is then given by 

<S = — E < 2 S [f(x,y)-f (x,y)] 2 I (4-2) 

s N 2 (x y J 

or 

5 = — E E E{f 2 (x, y)}- EE E{f(x, y) f(x, y)] 

3 n x y n x y 

+ — E E E{ f (x, y)} (4-3) 

n x y 

The first term above is the spatial domain autocorrelation 
function R(0, 0, 0, 0). The other terms may be evaluated by 
substitution of the reverse transforms yielding' 
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= R(0, 0,0,0) - — E D E 


/ r n2_i i 

jLw^O F ( w ) *b(x, y,w) 

£ E F(w ! ) b(x,y,w')J | +~2 E F ( w ) b (x, Y> w )j 

r m-i V 

E F(w') bCx^jW 1 ) 

L w'=0 J 


F(w') Mx^w 1 ) ; (4-4) 


Expanding the series and changing the order of summation gives 


I N -1 M-l 

$ = R(0,0,0, 0) — — T e \ E F(w)F(w J ) E E b(x, y, w) b(x, y, w') 


N 


M-l M-l 


.2 ] w= 0 w » = o 

w 4 w 1 


x y 


+ N E |w=0 ^ ^F(w)F(w')EE b(x,y,w) b(x,y,w') > (4-5) 


By the orthogonality of the kernel b(x, y, w) 


2 


N -1 M-l 


S = R(0, 0, 0, 0) - — ~ E \ E_ E F(w) F(w') 6 (w-w ') 
S n w- u w i = o 


I M-l M-l 

+ — 5 E < E S F(w) F(w') S (w-w 1 ) 


N 


w=0 w 0 


Thus, 


M-l 


<S =R(0, 0,0,0) - -^5 S E{F 2 (w)} 
S N w=0 


(4-6] 


(4-7a) 


or 
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1 M “ 1 2 

$ = C(0, 0,0,0) - ■ 2 c (w) (4- 7b) 

s w= 0 

For a given number, M, of samples to be' included in the transformation, 
the mean square error will be small if the variance of each of 
the transform samples is large. The transform which minimizes 
the mean square error is the Karhunen-Loeve transform of the 
first kind in which the eigenvectors are arranged in correspondence 
with the eigenvalues in descending order [30, 31] , 

For an image transform of the second kind, with zonal filtering, 
the reconstructed image is given by 

A N-l N-l 

f(x,y) = S 2 F(u, v) b(x,y,u, v) (4-8) 

u, v e M(u, v) 

where the transform domain indices are members of a set 
determined by a mask function M(u, v). By an analysis similar 
to that for an image transform of the first kind, it is found that 
the mean square error is of the form 

N-l N-l 

S s = C(0,0,0,0) --n- S 0 S 0 a (u, v) (4-9) 

u, v e M(u, v) 
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The transform which minimizes the mean square error is 
the Karhunen-L,oeve transform of the second kind for which the 
mask function corresponds, to the index pairs u, v that have the 
largest eigenvalues. 

Karhunen-Loeve Transform 

For the Karhunen-Loeve transform of the first kind, the 
minimum mean square error becomes 

M-l 

<? s = G(0, 0, 0, 0) — (4-10) 

N 

where X (w) represents the eigenvalues of the covariance matrix 
of the image. The operational procedure for performing zonal 
sampling with the Karhunen-Toeve transform of the first kind is 
simply to compute and code only the first M components of the -transform 
which are subsequently to be used in the inverse transform. . 

The minimum mean square reconstruction error for zonal 
sampling of the Karhunen- Loeve transform of the second kind is 

N-l N-l 

*S = C(0 ’°> 0 > 0 > - -^2 lo^Vfu.v) (4-11) 

u, v e M(u, v) 
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It should be noted that the mask function is generally not a 
simple rectangle in the transform domain. Exhibit 4-1 shows 
the ordering of the eigenvalues for a Karhunen-Loeve transform 
of the second kind. In this example the image covariance function 
is separable, and the vertical and horizontal element correlation 
is the same. Therefore, the eigenvalues corresponding to each 
line and column of the image are identical. The eigenvalue products 
give the same eigenvalues that would be obtained for a Karhunen- 
Eoeve transform of the first kind. However, there is no. simple 
and general ordering between X (u, v) and \(w). Thus, the eigenvalue 
ordering must be determined experimentally for a given Karhunen- 
Loeve transform. Figure 4-1 shows 16 by 16 element sampling 
masks for a 4:1 sample reduction for two values of the image 
covariance function. 

If a rectangular mask function 

u, v e M(u, v) if u < u^ ; v < v q rectangular mask (4-12) 

is employed for ease of implementation, the performance of the 
operation will not be optimum, but the degradation will usually not 
be too serious. Two other simple mask functions that could be 
employed are listed below: 
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w 

1 

_2 

3 

4 

5 

6 

n 

8 

(u,v) 

1,1 

1,2 

2,1 

1,3 

D 

D 

m 

2,2 

A (w) 

9.630 

1.735 

1.735 

.648 

.648 

.400 

.400 

.313 


w 

9 , 

10 

11 


13 

14 

D 

16 

(u,v) 

3,2 *. 

- 2,3. 

•4,2 



4,3 

| 


AM 

.117 

.117 

.072 

.072 

.044 

.027 

.027 

.017 


b. Ordering of eigenvalue products 


Exhibit 4-1. Karhunen-Loeve Transform Zonal Sampling Mask Generation 
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YC 
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Image covariance function 
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(XC) 


“ X, 


Figure 4-1 Karhunen-Loeve Zonal Sampling Masks 





(4-13) 
(4- 14) 

The hyperbolic mask most closely resembles the optimum 
mask determined by ordering the eigenvalues of the image 
covariance function. 

Fourier Transform 

Zonal sampling with the Fourier transform consists of 
sampling the lowest spatial frequencies in the transform domain. 

For the Fourier transform defined by Equation (2-20), the lowest 
spatial frequencies He in the four zones shown below: 




- 53 - 




Hadamard Transform 


With an ordered Hadamard transform, zonal sampling 
consists of sampling the transform domain samples with the 
low sequencies. These samples lie within a circular quadrant 
about the-origin in the transform domain. 

4, 2 Generalized Threshold Sampling 

Zonal sampling in the transform domain will provide small 
mean square error reconstructions of good subjective quality 
if the actual magnitude of a transform domain sample does not 
differ greatly from the standard deviation 0 (w) or <J(u, v). The 
difficulty with zonal sampling is that in most natural images there 
are many high spatial frequency samples lying outside the sampling 
zone that are of significant magnitude. In threshold sampling 
rather than determining a priori which transform domain samples 
are to be coded, the selection is made after the transform has been 
taken on a particular image. A threshold level is established 
a priori, or perhaps adaptively, and only those samples whose 
magnitudes are greater than the threshold are coded. If the 
threshold level is chosen a priori, based upon the probability density 
of the transform samples, the actual sample reduction factor for 
a particular image will be variable. As an alternative procedure 
the threshold level could be chosen so that a given number of 
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transform domain samples would be coded for a parti cular 
image. 

If transform domain threshold coding is -to be employed, 
the major question of interest is: What is the optimum image 
transform? In general the best transform is the transform which 
maps the image energy into the fewest transform domain samples. 
For a checkerboard image of half black and half white elements 
in each direction, the Hadamard transform is a very efficient 
transform since the image can be represented by only two 
transform domain samples. For natural images, the image can 
only be defined statistically, not deterministically. In such 
instances the optimum { minimum mean square error) transform 
is the transform for which the smallest number of samples have 
the largest variances. As stated previously, for a given class of 
images,- this transform is the Karhunen-Loeve transform. Thus, 
it is expected that the Karhunen-Loeve transform would exhibit 
the best minimum mean square' error performance for threshold 
coding for a given class of images. However, for a particular 
image of the class, another transform could provide better 
performance. 
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4. 3 Image Block Size Considerations 


For either zonal or threshold sampling in the transform 
domain, consideration must be given to the size of the image' 
block. From the standpoint of image energy it is best to make the 
image block size as large as possible in order to derive benefit 
from all element- to- element correlations within the image. 

For natural images, however, the correlation between elements 
separated by over 10 to 20 elements is usually relatively small- 
(see Appendix A). Therefore, little is lost in taking the image 
transform over smaller size blocks. This point is illustrated 
by Figure 4-2 which contains a plot of the percentage of tr an sform 
domain energy contained in the lowest one -fourth of the transform 
domain samples for a one dimensional Karhunen-Loeve transform 
as a function of block size. In this example the image covariance 
function is modeled as a Gauss -Markov process dependent only 
upon the adjacent element correlation factor, XC. As indicated 
in Figure 4-2 about 90% of the image energy is contained in the 
lowest one-fourth of the Karhunen-LiOeve transform coefficients 
for a block size of 16 by 16 elements. The percentage of energy 
contained in the low pass zone increases rather slowly for larger 
size blocks and decreases much more rapidly for smaller size 
blocks. It appears that a block size of about 16 by 16 elements 
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Percentage of Image Energy in Low Pass Zone, 


N/4 



4 8 


Figure 4-2 Effect of Image Block Size 
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Is a good compromise between maximizing the amount of image 
compression possible and simplifying the implementation of the 
transform coder* 

4*4 Comparison of Image Transforms 

A series of experiments has been conducted to determine 
the image coding performance of the Fourier, Hadamard, and 
Karhunen- Loeve transforms for natural images* As a result of 
the computational requirements of the Karhunen -Loeve transform, the 
image block size was limited to 16 by 16 elements* * 

Figure 4-3, 4-4, and 4-5 contain displays of the Fourier, 
Hadamard, and Karhunen- Loeve transforms* It should be noted 
that there is no apparent grid structure in the reconstructed 
image despite the block processing* 

Figures 4-6 to 4-10 illustrate the effect of zonal low pass 
filtering for the three transforms. In Figures 4-6 and 4-8 the filter 
pass band for the Fourier and Hadamard transforms is a circular 
zone in the transform domain* For sample reduction factors 
greater than 4:1 the 16 by 16 element grid structure becomes 
apparent because many of the transform samples that correspond to 


Examples of Fourier and Hadamard transform coding in larger 
size blocks are presented in Section 6* 
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b. Inverse transform of 
transform 


a. Threshold display of 
Surveyor box transform 



c. Threshold display of 
Girl transform 


d* Inverse transform 
of transform 



Figure 4-3 Fourier Transforms in 16 x IS Element Blocks 
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rirf * 



a* Threshold display of 
Surveyor box: transform 



b. Inverse transform of 
transform 


N 0^ 






c. Threshold display of 
Girl transform 


d, Inverse transform 
of transform 


figure 4-4 Hadamard Transforms in 16 x 16 Element Blocks 
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a, Threshold display of 
Surveyor box transform 


b. Inverse transform 
of transform 




c. Threshold display of 
Girl transform 


d* Inverse transform 
of transform 


not reproducible 


Figure 4-5 Karhunen-Loeve Transforms in 16 x 16 Element Blocks 


- 61 - 






a* 2:1 sample reduction 



b. 4:1 sample reduction 



Figure 4-6 Fourier Transform Zonal Sampling in 16 x 16 Element Blocks 
-- Circular Zone 
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b, 4:1 sample reduction 



c. 8:1 sample reduction 



Figure 4-7 Fourier Transform Zonal Sampling in 1 6 x 16 Element Blocks 
— Hyperbolic Zone 
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c. 8:1 sample reduction 



d. 12:1 sample reduction 


Figure 4-8 Hadamard Transform Zonal Sampling in 16 x 16 Element Blocks 

-- Circular Zone 
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c. 8:1 sample reduction d* 12:1 sample reduction 


figure 4-9 Hadamard Transform Zonal Sampling in 16 x 16 Element Blocks 
-- Hyperbolic Zone 
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Figure 4^igKarhunen-Loeve Transform Zonal Sampling in 16x16 Element Blocks 
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brightness changes in periods of 16 elements are excluded 
from the circular pass band* To prevent this grid effect a hyper- 
bolic shaped zone similar to the sampling mask of Figure 4- la was 
employed* The results shown in Figures 4-7 and 4-9 show a 
definite improvement in the elimination of the grid structure as 
compared to Figures 4-6 and 4-8. In all four images there is an 
expected loss of resolution. Figure 4-10 illustrates zonal low 
pass filtering with the Karhunen- Loeve transform. The filter is 
a mask passing those transform domain samples corresponding 
to the largest eigenvalues of the image covariance function. The 
reconstructed images do not show the 16 by 16 element grid 
structure p but there is some loss in resolution. Summarizing 
these results: for a given sample reduction factor the Karhunen- 
Loeve transform results in the smallest mean square error and 
the least image degradation from a subjective viewpoint* With 
a hyperbolic shaped pass band the Fourier transform is somewhat 
better than the Hadamard transform for both measures of image quality* 
Figures 4-11, 4-12, and 4-13 show the effects of threshold 
coding in the transform domain for the three types of image 
transforms. The quality rating of the three transforms from the 
standpoint of subjective quality is : Karhunen -Loeve, first: 

Hadamard, second; and Fourier, third. It should be noted that 
the sample reduction factors obtained for equivalent image quality 
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a, 5:1 sample reduction b, 10:1 sample reduction 



c* 20:1 sample reduction 


Figure 4-11 Fourier Transform Threshold Sampling in 16 x 16 Element Blocks 
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Figure 4“12 Hadamard Transform Threshold Sampling in 16 x 16 Element Blocks 
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c. 20:1 sample reduction 


Figure 4-13 Karhunen-Loeve Transform Threshold Sampling 
in 16 x 16 Element Blocks 
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are much higher for threshold sampling than zonal sampling. 

For a fair comparison, however, to account for position code 
bits, the sample reduction factors for threshold coding should 
be multiplied by a factor of about 0, 6 to 0* 8 to obtain 
equivalent bandwidth reduction factors. 

In summary the following conclusions can be drawn from 
these experiments: 

a. For both zonal and threshold sampling in the transform domain, 
the best transform is the Karhunen-Loeve transform, followed 
by the Hadamard transform, followed by the Fourier transform. 

b. Threshold sampling provides higher sample and bandwidth 
reduction factors than zonal sampling for all three transforms, 

c. The effect of a limited block size does not appear to be a 
serious problem either from the standpoint of image quality 
or performance. 

While the Karhunen- Loeve transform does appear to provide better 
performance than the Fourier and Hadamard transforms, the 
margin of performance is not too large. In view of the considerably 
greater amount of computation involved with the Karhunen-Loeve 
transform as compared to the Fourier and Hadamard transforms, 
its utilization will probably be limited. The following sections are 
therefore restricted in scope to the Fourier and Hadamard transforms. 
These sections present an analysis of transform domain quantization, 
a further discussion of image coding for bandwidth reduction, and a 
study of the error tolerance properties of image transforms. 
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5. Fourier and Hadamard Image Transform Quantization 

The dynamic range of Fourier and Hadamard transform 

domain samples in integer arithmetic is 1 to N A where N 

denotes the number of elements per line of the image and A is the 

maximum integer value of the amplitude of an image sample. 

If each transform domain sample were simply coded in a binary 
2 

code, log^(N A) bits would be required for each code word. 

For a 256 X 256 element image of 64 grey levels, each code wore* 
would be 22 bits in length. Even with threshold coding it would 
be unlikely that a significant bandwidth compression could be 
achieved for such large length code words. In order to achieve 
a bandwidth compression with transform coding it is necessary 
o recode, or quantize, each transform domain sample so that 
it may be represented by relatively short length code words. 

There are two basic approaches to this process: each sample could 
be quantized to the same number of levels, with the quantization 
levels possibly chosen according to a nonlinear scale; or the 
number of levels could be permitted to vary from sample to sample 
with a linear spacing of quantization levels. The latter approach 
will result in the most efficient coding, but the code words will be 
of variable length. This creates problems in data synchronization 
and channel coding for error correction. The former method can 
be adapted for relatively efficient constant word length coding. 
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In the subsequent discussion only the first method of quantization 
is considered. 

5. 1 Quantization Scales 

In the quantization process let the transform sample 

component (amplitude, real part, imaginary part, magnitude, 

or phase) to be quantized be represented by the function 

F (u, v). The range of the sample component is assumed to be 
c 

broken up into K positive and K negative bands separated by 

th 

quantization levels Q^(j = 0, ±1, ±2, •••, ±K). The zero 
quantization level and the upper and lower quantization levels 
are assigned the values 

(5- la) 
(5- lb) 
(5-lc) 


Q 0 =0 

^ NA 
°K = ~2“ 


Q 


~K 


-NA 

2 


where A represents the maximum value of a sample of the original 
image of N by N -elements. If a transform component falls in a band 
bounded by quantization levels Q^ ^ and Q^, the component is 
quantized, and subsequently reconstructed, to the value F^(u, v) 
which lies within the band. The relationship between quantization 
and reconstruction levels is given in Figure 5-1. 

Quantization and reconstruction levels are logically chosen 
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reconstruction 

levels 
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NA 

2 



F 


-2 


A 


Q 


-K+l 


T 


-i 

A_ 


Q 


0 



K-l 




K 

-A- 


K-l 


NA 


quantization 

levels 


Figure 5-1 Quantization and Reconstruction Levels 



to minimize the effects of the quantization error introduced by 
the amplitude truncation of samples. Table 5-1 lists some error 
criteria for the selection of quantization and reconstruction 
levels. The quantization error criterion depends upon the appli- 
cation of a reconstructed image. The principal consideration is 
whether the image is to be used for subjective viewing or photo- 
metric measurements. 

For subjective viewing the relative spatial error criterion 
listed in Table 5-1 provides 'an indication of image quality. This 
relative spatial error criterion is predicated upon the fact that 
incremental brightness changes in the reconstructed image are much 
more noticeable if the brightness level is low than if it is high. 

Thus, to minimize the relative spatial error, the density of 
quantization levels in the spatial domain should be greater at the 
lower amplitude levels. But, since the brightness of every point 
of a reconstructed image is a function of the amplitude of a single 
transform sample, then by the same reasoning, the density of 
quantization levels should be greater for low level transform samples. 
From psychophysical tests, it is known that the human viewer is 
very sensitive to the location of high frequency brightness transitions, 
but relatively insensitive to their actual magnitude. In fact 
images which have been "crispened" by high pass filtering often 
appear, preferable to the original image. From this characteristic 
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TABLE 5-1 


Quantization Error Criteria 


Cumulative mean -square 
spatial error 

N-l N-l . 

2 2 yM(x,y)] 2 

X=0 y=0 L ' J 

Cumulative mean square 
transform error 

N-l N-l 

2 2 [ F (uv v)-F(u, v)l 

u=0 v=0 L J 

Cumulative spatial error 

N-l N-l 

2 2 f(x,y)-f(x, y)| 

x=0 y=0 

Cumulative transform error 

■ N-l N-.l ^ ■ , 

2 2 v >] 
u=0 v=0 

Relative spatial error 

if(x, y)-f(x-, y)i 

|f(x. y) | 

Relative transform error* 

!f(u, v)-F(u, v)l 
|F(u,v)| 


F(u f v) = quantized value of F(u, v) 


?(x. y) 


inverse transform of quantized value of F(u,v) 
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of subjective viewing it would seem that the density of quantization 
levels, at low transform sample amplitudes, should be greater 
at the higher spatial frequencies than at the lower spatial 
frequencies. 

If photometric measurements are to be made on an image 
the cumulative mean square spatial error is a common fidelity 
criterion. For a mean square error criterion the quantization 
levels in the transform domain must be selected to minimize 
the cumulative mean square error in the spatial domain. Let 

N-l N-l f _ \ 

6 S " ^2 xio S 0 E | [f(x * y) - f(x ’ y)] f (5 ‘ 2) 

represent the cumulative mean square spatial domain error 
where f (x, y) is the image reconstruction from the quantized 
transform samples, F(u, v). For a Fourier or Hadamard transform 

N-l N-l 

£ (x, y) = Ij_ q F(u,v) b(x,.y,u,v) (5-3) 


and 


Nr 1 N-l ^ 

f(x, y) = Tl 0 Tj 0 F(u, v) b(x, y, u, v) 


(5-4) 


Hence, the spatial domain mean square error can be written as 
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N-l N-l 

6 = T, S„ E< 

s * x= 0 y= 0 


N-l ' N-l 

5=0 v S Q [F.(u,v)-F(u,v)] 


•b(x, y,u,v) 


Expanding the integrand yields 


N-l N-l 

6 = — o EL E„ e< 


N 


,2 x= 0 y= 0 


N-l N-l N-l N-l 

u 5 0 V S 0 E) E [F(u,v)-F(u, v)3 

u' = 0 v*=0 


■ [F(u' } v'J-FJu 1 , v 1 )] b(x, y, u, v) b^y^jv') 


Rearranging the order of summation gives' 
N-l N-l N-l N-l 

E 

u 1 = 0 v f — 0 


w-l iN-l W-l 1 

= — i iMovio S S E { \LF(u,v)-F(u,v)] 
N u ( = 0 v ! = 0 1 


N-l N-l 

[F(u , ) v , )-F(u , ,v , )]‘ Ij q E q b(x, y,u, v) b(x, y, u’, v 1 )] 


As a result of the orthogonality of the transform kernel, 


6 - 2 E ( [F(u,v)-F(u, v)] [F(u , ,v i )-F(u , ,v 1 )] 

N u 1 v l i 

• 6{u-u ! , v-v 1 ) ! 


6 q = "1 JS E \ [F(u, v)-F(u, v)] 2 

s u v I 


(5-5) 


(5-6) 


(5-7) 


(5- 8a) 


(5- 8b) 
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The cumulative mean square error in the spatial domain of 
the reconstructed image is therefore equal to the cumulative 
mean square error in the transform domain. Minimization of 

can then be accomplished by minimization of the mean square 
error 


<2 (u, v) 


E‘ < [F(u, v)~F(u, v)]* 


(5-9) 


in the transform domain for all spatial frequencies. For the 
Hadamard transform, quantization and reconstruction levels 
for the transform sample amplitude must be found to minimize 
<$ (u, v). In the case of the Fourier transform the mean square 
error of the real and imaginary, or phase and magnitude , 
components of a transform must each be minimized. The mean 
square error of a transform component may be written in explicit 
form as 


& (u, v) = <S + (u, v) + <$ (u; v) 

in which 


(5-10) 


TT Q-(U,v) 

J 2 

$ (u,v) = £ J [F (u, v)-F.(u,v)] p(F ) dF (5-lla) 

J Qj.iK-) c J c ° 

and 
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(5- lib) 


-K Q j (U ’ v) 

5 (u,v) = S J* ’ [F (u,v)-F.(u,v)] p(F ) d F 
j=-l Q. + 1 (u,v) C 3 c c 

where p(F^) is the probability density of the transform sample 
component to be quantized. If p(F^) is a symmetrical probability 
density about Q Q = 0, then <S + (u, v) equals <$ (u,v), and the 
quantization rule determined by the minimization of <S + (u,v) 
is the same as that determined from <$ (u,v). 

The optimum placement of the quantization and reconstruction 
levels to minimize the mean square error of a quantized signal 
has been determined by Panter and Dite [40 ] . The reconstruction 
levels should be located midway between each pair of quantization 
levels , Thus 


Q.(u, v) -f Q (u, v) 
Fj(u, v) = -J 


(5-12) 


The quantization levels can be determined to a good approximation 
[40] from 


. NA 


NA 


Q.(u,v) = 


PfF c ] 


NA 

2 


J p{F ] d F 


d F 


(5-13) 


- 80 - 



Three cases of interest for quantization of the Fourier and 


Hadamard transforms are listed below. 


Uniform distribution: 

P [F c 3 = NA 
Q j( u ' v > = ■> 2K 


Rayleigh distribution:, 

F 

p{F c J = 9 exp { - 


2 , 

a (u, v) 


2 , 

F c (u,v) 
20‘ 2 {u,v) 


Q jKv) 


J F (u,v) ] 

f [F (u,v)] exp \ — } d F 

0 c 


NA 

NA r J 2K 
2 


60 2 (u,v) 


NA i 

.”F” " 3 

J [F (u,v)] exp 


F 2 (u, v) 
6 0 2 (u, v) 


d F 


Gaussian distribution: 


p{ F ] = [2tt0 (u,v)] exp \ - — ~ 


2 , 

F (u,v) 


2CT 2 (u, v) 


NA 


.NA 
• 3 2K 


exp 


Q j (U ’ V)= "~NA- 


F 2 (u, v) 
c 

2 

6 0 (u,'v) 


d F 


exp 


F 2 (n, V ) 

> d F 

60 (u,v)J 


(5- 14a) 
(5- 14b) 


(5-15a) 


’ (5- 15b) 


(5-1 6a) 


(5- 16b) 
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The quantization scales determined by Equations (5- 15b) 

and (5-1 6b) for the Rayleigh and. Gaussian distributions have 

the desired subjective property that the quantization levels are 

more closely spaced at the lower quantization levels, and a _re 

more closely spaced at the higher spatial frequencies for which the 
2 

variance 0 (u, v) is smaller. Unfortunately, the quantization levels 
are nonlinearly related to the sample variance. Hence, it 
becomes necessary to compute a separate quantization scale for 
each transform sample. 

There are two other scaling laws- -the Gaussian error 

function and the logarithmic- -that have the same general 

characteristics as the optimum mean square error quantizer, but 

for which the quantization levels are linearly related to the 

sample variance, ha the Gaussian error function quantizer the 

quantization levels are selected so that when the probability density 

Z 

of transform samples is Gaussian with variance Q .(u, v), the 
probability that a transform sample is quantized to a given 
reconstruction level is the same for all levels. This results in 
a uniform entropy for all reconstruction levels, and therefore a 
constant word length code may be used for each quantized sample. 
The quantization levels are given by the solution of the equation 
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2 

K 


Q. , 

J J 7 ■■ 1 

Qj_ 1 V 2rr a z (u, v ) 


exp 


J 


NA 

2 


0 V 2tt a 2 ( Uj v ) 


exp 


F^(u, v) 

20 2 (u, v) 

— 2 for j= 1,2, • • • ,K-1 

F (u,v) | 

S — > dF c 

20 (u,v) 



(5-17) 


For “ — large the denominator approaches one -half and the 
scaling law can be expressed in terms of the Gaussian error function 
as 


1 

2K 


1 > erf 


Q. 


Jz O (u, v) 


— erf 


AM 

JT 0 (u, v) 


(5-18) 


where 


erf{x} s* — J exp{-Z 2 } dz 
- 0 


The logarithmic quantizer- obeys the function 


An 



q: 

W(u, v) 


K 


An 



, for j = 0, 1, 2, 


,’K-l 


(5-19) 


in the positive quadrat and the inverted and reversed function 
in the negative quadrant where W(u, v) is a spatial frequency 
weighting function. The quantization levels are approximately 


- 83 - 



Q. = 

J 


W(u, v) 


for j = 0,1,2,°“*, K-l 


(5-20a) 



(5-20b) 


A convenient implementation of the logarithmic quantizer can 
be realized by adding plus one to each sample component and then 
taking the logarithm. The resulting continuous function can 
then be quantized linearly. 

Figure 5r-2 shows the relationship between the quantization 

levels set by the optimum, Gaussian error function, and logarithmic 

quantizers when the probability density of the transform domain 

2 

samples is Gaussian with variance O' (u, v). This figure indicates 
that Gaussian error function scale is a reasonably good approximation 
to the optimum scale for a transform sample maximum standard 
deviation in the range of about 1, 500 to 4, 000. 

5, 2 Quantization Experiments 

A series of experiments has been conducted to assess the 
effects of quantization of Fourier and Hadamard transform domain 
samples. In these experiments the transform domains were 
quantized and reconstructions were obtained of the quantized samples. 
The cumulative root mean square quantization error 
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LEVEL NUMBER 


Figure 5-2 Comparison of Optimum and Gaussian Error Function 
Quantization Scales 
-85- 




6 

Q 


1 


(5-21) 


2 5 ? 
1 y 


It 


N 


J 


was measured for each quantized image. In addition the 


difference function 


d(x,y) = |f(x, y}-f (x, y) | 


(5-22) 


was formed to indicate the spatial correlation of errors. 

Figures 5-3 and 5-6 show the effects of quantization on 
the Fourier and Hadamard transforms, respectively, for a 
Gaussian error function quantizer with 64 quantization levels. 

In these experiments the transform domain variance function was 
modeled as 


where S and p are the maximum and spread variance parameters. 
A computer search procedure was developed to determine the 
best values of S and p to minimize the quantization error S _ „ 

The reconstructions in Figures 5-3 and 5-6 were made on images 
quantized with the values of S and p giving a minimum value of 
Figures 5-4, 5-5, 5-7, and 5-8 illustrate the effect of an 
incorrect choice of the variance parameters S and p. There is a 
broad range in the values of S and p which provide good quality 
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amp . par, = 91 ,000 
spread par, = 4,000 



b, difference 
R*M,S, error = 2*2 


Surveyor box 



c, amp. par, = 174,000 
spread par, = 1,500 



R,M.S. error =2,9 


Surveyor boom 



Figure 5-3 Fourier transform quantization: examples of correct parameter 
scaling — Gaussian error function quantizer , 64 levels 
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a, amp. par. = 350,000 
spread par, = 4 , 000 


b. difference 
R. M.S. error - 11,9 



too Large amplitude parameter 



too small amplitude parameter 


Proper spread parameter 


Figure S-4 Fourier transform quantization: examples of incorrect amplitude 

parameter scaling — Gaussian error function quantizer, 64 Levels 
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too large spread parameter 



c, amp, par, = 91000 
spread par* - 1000 



d t difference 
R . M . S . error - 5,3 


too small spread parameter 


Proper amplitude parameter 


Figure 5-5 Fourier transform quantization' examples of incorrect spread 
parameter — Gaussian error function quantizer, 64 levels 


- 89 - 




Surveyor box 







e, amp, par* = 400 
spread par* = 10,000 

Surveyor footpad 


f * difference 
R* M * S * error = 1.3 


Figure 5-6 Hadamard transform quantization: examples of correct parameter 
scaling -- Gaussian error function 
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too large amplitude parameter 



c* amp. par* - 200 

spread par* = 11,100 



too small amplitude parameter 


Proper spread amplitude 


Figure 5-7 Hadamard transform quantization: examples of incorrect amplitude 
parameter scaling - Gaussian error function quantizer, 64 levels 
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too large spread parameter 




c. amp* par* - 475 
spread par* - 4 000 


d* difference 
R * M * S ■ errors 4*2 


too small spread parameter 


Proper amplitude parameter 


Figure 5-8 Hadamard transform quantisations examples of incorrect spread 
parameter -- Gaussian error function quantizer, 84 levels 
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reconstructions . Experiment jb previously reported [27] show 
the effect of using 32 and 16 quantization levels picked according 
to the Gaussian error function quantization scale. The quantization 
error is noticeable for 32 levels and quite bad for 16 levels. 

In summary of the quantization experiments, it has been 
found that good quality Fourier and Hadamard transform 
reconstructions are possible when the transform samples have been 
quantized to as few as 64 levels using the Gaussian error function 
quantization scale. 
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6. Fourier and Hadamard Image Transform Bandwidth Reduction 
A sample reduction, and a subsequent bandwidth reduction, 
by proper coding, are possible by coding the Fourier or Hadamard 
transform of an image rather than the image itself. This sample 
reduction is obtainable because, as a result of the element- to* 
element correlation in the image, many of the transform domain 
samples are of extremely low magnitude and may be deleted 
from the image reconstruction without seriously degrading the 
quality of the reconstructed image* 

The process of selecting samples for inclusion in the image 
reconstruction can be conveniently analyzed from the viewpoint 
of two dimensional sampling. Figure 6-1 illustrates a generalized 
block diagram of a transform sampling system. The forward 
transform of an image, F(u,v), is multiplied by a two dimensional 
sampling function, S(u,v), which takes on the values zero or one 
according to some a priori or adaptive rule* The sampled image 
transform, F s (u,v), is 

F g (u,v) = F(u,v) S(u,v) (6-1) 

The reconstructed image, f^(x,y), is then the reverse transform 

of F (u,v), Thus, 
s 

N-I N- 1 

f s (x, y) - n S 0 F(u, v) S(u, v) b(x,y, u, v) (6-2) 
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F( u , v ) F (u. 


original 

image 



FORWARD 

TRANSFORM 


f(x,y) 


o 

S (u , V ) 

sampling 

function 


REVERSE 

TRANSFORM 


Figure 6-1 Transform Domain Sampling 


In the case of the Fourier transform, as a consequence of 
the frequency translation theorem, the reconstructed image 
can be expressed as a spatial convolution, denoted by the 
symbol © , of the original image and the inverse Fourier 

transform, s(x,y), of S(u,v). Thus, 

f(x, y) = f fl (x, y) ® s(x, y) (6-3) 

It should be noted that this result does not hold for the Hadamard 
transform since the Hadamard transform does not possess a 
sequency translation property. 

Table 6-1 lists three basic transform sampling methods. 

With the random sampling method the sampling function, S(u, v), 
assumes the value 0 or 1 according to some probability 
distribution p(u, v) over the transform domain. Experiments have 
been performed in which one-half of the Fourier transform samples 
have been randomly discarded independent of their location in the 
transform domain. The resultant reconstructions were of poor 
quality due to errors in deleting large magnitude low frequency 
samples. Several variations were attempted in which more of the lower 
spatial frequencies were included, but die results were not 
particularly encouraging. It appears that at most a 2:1 sample 
reduction can be obtained by random sampling at the cost of 
moderate degradation. 
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TABLE 6-1 


Classification of Transform Sampling Methods 


Description 

Sampling 

Function 

S(u,v) 

Conditions 

Random 

sampling 

1 

with probability p(u F v) 

0 

with probability 1 - p(u , v) 

Zonal 

sampling 

1 

u,v in sampling region 

0 

u,v not in sampling region 

Threshold 

sampling 

1 

j F(u,v) j > M t (u,v) 

0 

J F(u,v) [ < M^{u , v) 
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The zonal and threshold sampling techniques aTe 
discussed in the following sections. 

6. 1 Zonal Transform Sampling 

In most scenes of interest the energy in the Fourier 
transform domain tends to be clustered toward the lowest spatial 
frequencies. Similarly, the Hadamard transform domain energy 
is greatest at the low sequencies* For example, in the three 
Surveyor spacecraft scenes, 95% of the image energy in the 
Fourier transform is contained in 1% or less of the transform 
samples [27]. 

With an image energy distribution clustered at the low 
spatial frequencies or sequencies, the most obvious means of 
conserving bandwidth is simply to not transmit the high spatial 
frequency or sequency samples. Discarding the high spatial 
frequencies or sequencies is equivalent to passing the image 
through a circular, zonal, low pass filter; the result is a loss of 
focus. If some degree of resolution loss is acceptable, zonal low 
pass filtering of the transform domain yields relatively large 
bandwidth reductions. 

Zonal low pass filtering of a sequency ordered Hadamard 
transform is equivalent to multiplying the transform samples by 
the sampling function 
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Hadamard transform: 


2 2 2 

S(u, v) = 0 if u + v > R 
S(u,v) s 1 otherwise 


u,v - 0, l, * * ' , N-I 


For the Fourier transform the sampling function is 
Fourier transform: 


(6-4) 


S(u P v) = 1 
S(N-l-u, v) = 1 
S(u, N-l-v) = 1 
S(N- 1 -u, N-l-v) - 1 
S(u» v) = 0 


A 


Z 2 1 

if u + v > R 

^,v - 0, l p * - * , 


othe rwise 

J 


(6-5) 


Figure 6-2 shows the effects of Fourier and Hadamard transform 
zonal low pass sampling of the Surveyor box scene over the full 
frame of 256 by 256 elements* These experiments support the 
widely known fact that the high frequency and sequency brightness 
transitions are important even though they are relatively few in 
number and contain a low proportion of the image energy* The 
image degradation tends to be more noticeable for zonal filtering 
of the Hadamard transform than the Fourier transform for the 
same sample reduction factor because of the rectangular shape of 
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a, 4:1 Fourier b. 4:1 Hadamard 

R = 143 



c. 8:1 Fourier d. 8:1 Hadamard 

R = 101 



e, 16:1 Fourier f, 16:1 Hadamard 

R = 71 


Figure 6-2* Circular Low Pass Zonal Fourier and Hadamard Transform 
Sampling of Surveyor Box over Full Frame of 256 x 256 
Elements, unquantized transform. 
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the two dimensional Hadamard reconstruction waveforms. The 
eye is very sensitive to the presence of sharp brightness transitions 
within an image. With the Hadamard transform all transitions 
occur within one element, whereas in the Fourier transform the 
brightness transitions are spread over many elements since the 
reconstruction waveforms are two dimensional sinusoids* 

If the zonal low pass filter has square rather than circular 
boundaries in the transform domain, it is possible to produce a 
low pass version of the original by the simple expedient of spatial 
averaging of elements in the original image* In this case the 
complexities of the transform operation would probably not be 
warranted if a low pass reconstruction is acceptable. 

It has been conjectured that to produce a subjectively pleasing 
image, the eye only requires the low spatial frequencies of an 
image signal to provide the overall grey scale and the high spatial 
frequencies to provide the edge transitions* The mid-spatial 
frequencies are assumed to play a minor part in the reconstruction 
of an image* This conjecture has been tested by sampling the 
Fourier and Hadamard transforms of an image with a circular 
zonal rejection filter with the characteristic functions 
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Fourier transform: 


S(u,v> = 0 


S(N-l-u.v) = 0 



S(u, N- 1 -v) - 0 


S(N- 1 -Uj N-l-v) = 0 


( 6 - 6 ) 


S(uj.v)' " 1 otherwise 


Hadamard transform: 



(6-7) 


S(u f v) = 1 Otherwise 

Figure 6-3 shows the effect of band rejection filtering on the Surveyor 
box scene. The image quality appears to be somewhat degraded as 
compared to the results of Figure 6-2 with a simple circular zonal 
low pass selection of transform samples. 

In the development of zonal transform sampling techniques 
presented in Section 4 for transforms in 16 by 16 element blocks, it 
was found that a hyperbolic zone suppressed the grid effect for high 
sample reduction factors better than a circular zone. The effect of 
the use of a hyperbolic zone of full size, 256 by 256 point, images is 
shown in Figure 6-4. These reconstructions show a perceptible 
improvement to their counterparts in Figure 6-2 for a circular zone. 
For natural images it appears that a hyperbolic zone matches the 
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a* Fourier b, Badamard 

Rj =100; R 2 =233,8 

R 3 ® 256 



c« Fourier d* Hadamard 

^=125; R z o 245*5 

R 3 * 256 


4:1 sample reduction 


Figure 6-3 Circular Band Rejection Zonal Fourier and Hadamard 
Transform Sampling of Surveyor Box over Full Frame 
of 256 x 256 Elements, unquantized transform. 
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"si 4 *! 


Figure 6-4 Hyperbolic Low Pass Zonal Fourier and Hadamard Transform 
Sampling of Surveyor Box over Full Frame of 256 x 2S6 
Elements, unquantized transform. 


a, 4:1 Fourier 


c. 6:1 Fourier 


b. 4:1 Hadamard 


d* 3:1 Hadamard 


REPS' 
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energy distribution better in a spatial frequency or sequency ordered 
transform plane than does a circular zone. 

6. 2 Threshold Transform Sampling 

The difficulty with zonal transform sampling is that large 
magnitude transform samples often are included in the rejection zone 
and therefore deleted from the reconstruction. This problem can be 
overcome by the threshold sampling technique in which samples whose 
magnitudes are greater than a pre- specified threshold level are 
included in the image reconstruction independent of their position 
in the transform domain. 

Figure 6-5 and 6-6 are plots of the percentage of transform 
domain samples lying below a magnitude threshold level for the Fourier 
and Hadamard transforms. Maps showing the location of transform 
samples exceeding the threshold level for the Fourier and Hadamard 
transforms are shown in Figure 6-7. It should be noted that the large 
magnitude samples tend to be located at the lower spatial frequencies 
or sequencies. But many high spatial frequency and sequency samples 
exceed the threshold. In low pass zonal filtering these transform domain 
samples would not have been included in the image reconstruction. 

Figure 6-8 to 6-11 show the effects of threshold coding in the 
transform domain for the Fourier and Hadamard transforms. Each 
transform domain has been quantized to 64 levels per transform sample 
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0 4 8 12 16 20 25 30 35x100 

THRESHOLD LEVEL 

Figure 6-5 Number of Fourier Transform Samples below Threshold 
versus Threshold Level 

-inf;- 




a. Fourier b. Hadamard 

5:1 sample reduction 



c, Fourier d. Hadamard 

10: 1 sample reduction 



e. Fourier I, Hadamard 

20:1 sample reduction 

Figure 6-7 Maps of Fourier and Hadamard Transform Samples above 
Threshold for Surveyor Box over a Full Frame of 256 x 256 
Elements . 
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a. 5; 1 sample reduction 


b . Difference 
R.M ,S. error = 3,6 




c. 10:1 sample reduction 


d. Difference 
R . M , S , error = 3,8 



e, 20:1 sample reduction 


f. Difference 
R , IV! , S . error - 4,7 



Figure 6-8 Fourier Transform Threshold Coding: Effects of Thresh- 
olding for Surveyor Box over a Full Frame of 256 x 256 
Elements , quantized transform. 
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a. 10:1 sample reduction 


b . Difference 
R . M . S , error * 3,7 


Surveyor Boom 



Surveyor Footpad 


Figure 6-9 Fourier Transform Threshold Coding: Effects of Thresh- 
olding for Surveyor Footpad and Boom over a Full Frame 
of 256 x 256 Elements , quantised transform. 
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c. 10:1 sample reduction 


d. Difference 
R. M.S , error - 3,9 




e, 20:1 sample reduction 



f* Difference 
R, M * S- error = 4, S 


Figure 6-10 Hadamard Transform Threshold Coding: Effects of Thresh- 
olding for Surveyor Box over a Full Frame of 256 x 256 
Elements, quantized transform. 


- 111 - 





a, 10:1 sample reduction 


t , Difference 
R - M , S , error = 4.8 


Surveyor Boom 



Surveyor Footpad 


Figure 6-11 Hadamard Transform Threshold Coding: Effects of Thresh- 
olding for Surveyor Footpad and Boom over a Full Frame 
of 256 x 256 Elements, quantized transform* 
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component according to the Gaussian error function quantizer scale. 
Thus, the image reconstructions exhibit the joint effects of sample 
deletion and quantization. Difference pictures are displayed to illustrate 
the spatial distribution of errors. Also the cumulative average mean 
square error has been measured for each reconstruction. From these 
experiments it can be concluded that the Fourier and Hadamard 
transforms both provide good quality reconstructions for sample reduc- 
tion factors of 5:1, Some image degradation is noticeable for a sample 
reduction factor of 10:1. 

In order to achieve a bandwidth reduction for digital image 
transmission with transform domain threshold coding it is necessary 
to code the position of the samples exceeding the threshold level. 

There are a variety of ways of position coding that could be employed. 
The simplest conceptually would be to code the coordinates of each 
significant transform sample. Higher coding efficiency can be obtained, 
however, by coding the number of non- significant samples between 
significant samples. This scheme, called run length coding, has been 
used quite successfully in the spatial domain for black or white pictures* 
To achieve a short position code length, runs are usually restricted in 
length to some maximum value, normally a power of two. By including 
a line synchronization code group it becomes unnecessary to code the 
line number. Another advantage of the employment of a line synchro- 
nization code is that it prevents the propagation of channel errors over 
more than one line. 
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A run length coding procedure for Fourier and Hadamard 
transform threshold coding has been implemented on a general 
purpose digital computer. The coding procedure is "fail safe" in 
that every transform domain sample is coded for a zero level threshold 
there is no truncation of transform samples. The basic properties 
of the run length coding procedure are outlined below: 

a. The first sample along each line is coded regardless of its 
magnitude. A position code of all zero bits is affixed to the 
amplitude code to compromise the line synchronization code 
group. 

b. The amplitude of the second run length code word is the coded 
amplitude of the next significant sample. The position code is 
the binary count of the number of samples of the significant 
sample from the previous significant sample* 

c. If a significant sample is not encountered after scanning the 
maximum run length of samples, the position code bits are 
set to all ones to indicate a maximum run length. 

A simple code to implement this run length coding procedure is 
given below. 
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This run length coding procedure for transform threshold coding 
has been tested for the Surveyor box and boom scenes. As expected, 
the run length coding does not introduce any reconstruction errors* 
The effect of channel errors on position bits is considered in the 
next section. Table 6-2 shows the bandwidth reduction factors 
obtained for these test scenes as a function of the sample reduction 
factor. In all cases the run length code employed four position 
bits and runs were truncated in length to 14 samples. Better per- 
formance could, no doubt, be obtained if the number of position code 
bits were tailored to match the run statistics. 
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TABLE 6-2 


Experimental Bandwidth Reduction Achieved by Fourier 
and Hadamard Transform Threshold Coding 
for Surveyor Box Scene 



Fourier Transform 

Hadamard Transform 

Sample 

Reduction 

Number of position bits 

Number of position bits 

3 

4 

S 

6 

3 

4 

5 

6 

S : 1 

2. 6 

3,1 

3.2 

3.2 

2.2 

2.6 

2,6 

2.5 

10 : 1 

3,3 

4, a 

5.7 

6.0 

3.0 

4,3 

4,9 

S.l 

20 : 1 

3,9 

6.6 

8.7 

10.0 

3.3 

5.4 

6.9 

7,7 


7* Fourier and Hadamard Image Transform 
Channel Error Tolerance 

A major concern of communication system designers is the 
susceptibility of data to noise interference* It is important, then, to 
study the effects of noise on the image transform coding communication 
system* The inherent 1 'error averaging” property of transform coding 
combined with error correction coding of specific transform samples 
provides a means of image coding for which channel errors are less 
deleterious than for conventional spatial coding of an image* This 
property, of course, is predicated on the assumption that the 
particular transform used tends to compact image energy in a few 
number of coefficients in the transform domain* 

In most digital communication systems the code alphabet 
consists of two symbols which are subject to perturbations in the 
channel, and these perturbations introduce random noise at the 
receiver. The binary symmetric channel is used as the noise model 
in the study of channel effects on image transform coding* The 
classical representation of such a communication channel is given 
in Figure 7-1, where the probability of receiving an incorrect symbol 
is p regardless of which symbol is transmitted* 
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Figure 7-1 Model of a Binary Symmetric Channel 
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7. 1 Channel Noise Effects 


An intuitive justification for transmitting the transform 
rather than the spatial domain of an image is the fact that channel 
noise introduced in the transform of an image tends to be distributed 
over the entire reconstructed image. Consequently, the noise 
manifests itself as a combination of low order orthogonal functions 
in the image due to noise introduced in the large amplitude coefficients 
of the transform domain. If the Fourier transform is used, the noise 
presents itself as a low frequency effect and if the Hadamard trans- 
form is used, the effect is low sequency corresponding to non-periodic 
checkerboards of a low number of zero crossings. Finally, if the 
Karhunen^ljoeve transform could be used, the noise introduced in the 
large valued coefficients would correspond to those orthogonal func- 
tions representing the largest eigenvalues and matching the original 
image closest in a mean square error sense. In all cases, since the 
eye is more sensitive to the high frequency ,r salt and pepper" effect 
of channel noise in the spatial domain, the same channel error rate 
in the transform domain is somewhat less offensive. Figure 7 -2a 
shows a mid- grey scene after having passed through a channel with 
probability of error of Pe - 10 Figure 7- 2b is the Fourier trans- 
form of the output of the same channel whose input was the Fourier 
transform of the mid-grey scene. Figure 7-2c is the same experiment 
replacing the Fourier transform with the Hadamard transform* All 
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a. BSC noise in spatial domain b* Fourier transform of BSC 

noise in Fourier domain 



c, Hadamard transform of BSC 
noise in Hadamard domain 


Figure 7-2 Binary Symmetric Channel with Error Rate Pe = 10 
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three scenes have the same error rate but the induced noise energy 
is distributed quite differently. A quantizing and coding method can 
be developed to take advantage of the inherent high frequency or 
"salt and pepper if noise immunity that transform domain coding 
offers. As a first step in this direction a requirement will be made 
that each quantum level occur equally likely as any other quantum 
level. This quantization criterion will guarantee that each code word 
is equally likely to occur and will avoid any unexpected noise biasing, 
since the binary symmetric channel affects each code bit, and there- 
fore each code word, independently of all others. Such a quantization 
requirement results in the quantization rule employed in the earlier 
sections of this report. As was mentioned earlier, such a scheme is 
sub-optimum with respect to quantization noise error, but is better 
suited for channel noise immunity. 

A sequence of computer noise simulation experiments have 
been conducted in order to verify the concepts presented earlier. 

Figures 7-3, 7-4, and 7-5 present the results of the simulation where 
three different noise rates: Fe = 10 ^ Fe =10 , Pe = 10 r were 

introduced into the spatial, Fourier, and Hadamard domains respectively. 
In addition the difference pictures are included for visual purposes. 

The "salt and pepper" effect is quite evident in Figure 7-3 for spatial 

-4 

domain errors. For errors less than Pe = 10 the transform 
domains indicate little or no degradation while a few errors are still 
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Figure 7-3 Spatial Domain Coding Effects of Channel Errors 
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Figure 7-4 Fourier Transform Coding Effects of Channel Errors 
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Figure 7-5 Hadamard Transform Coding Effects of Channel Errors 
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evident in the spatial domain noise, However, for larger noise rates 
the low order orthogonal functions that make up the respective trans- 
formations tend to swamp out the reconstructed image. This can be 
explained by the fact that the absolute, as opposed to relative value 
of a bit error is much larger in the regions where the transform 
coefficients (eigenvalues for the Karhunen- Loeve transform) are large 
in the transform domain. This explains the effect in Figures 7-4e 
and 7-5e, Further demonstration of this effect was presented in 
reference [27] where it was shown that by protecting certain areas 
of the transform domain from noise effects, large improvements in 
noise immunity could be obtained. This suggests an error correction 
procedure, a simulation of which is presented in the following section. 
However, before developing some error correction techniques, it is 
instructive to investigate the effects of a noisy channel on thresholded 
transform domains in order that both bandwidth reduction and noise 
immunity be combined. Figures 7-6 and 7-7 present results of such 
a simulation in which a threshold has been selected to provide a 5:1 
sample reduction ratio. Again the difference pictures are presented 
for visual evaluation purposes. The noise effects now include run 
length errors in the transform domain which manifest themselves as a 

unique type of one dimensional blurring in the reconstructed images. 

_4 

Again, noise with errors less than Pe - 10 tend to be averaged out 
due to the reconstruction process. The threshold coding technique 
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b. Difference 



d , Difference 





Figure 7-6 Fourier Transform Threshold Coding 
Effects of Channel Errors S*R. =5:1 
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Figure 7-7 


Hadamard Transform Threshold Coding 
Effects of Channel Errors, S.R, =5:1 
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requires coding for position information. These code words should 
be uniformly distributed so that unexpected noise biasing does not 
occur in the position code as well as data code. More sophisticated 
coding techniques might be pursued in this area, 

7, 2 Error Correction Transform Coding 

As a result of the statistical regularity of samples in the 
transform domain, a smaller amount of error correction in this 
domain will yield a better noise immunity than the same amount of 
error correction in the spatial domain. The nature of the quantization 
law is such that errors in certain positions of the transform domain 
are much more bothersome than in other positions due to the large 
statistical variance of samples at these coefficients. Therefore, it 
is natural to develop an error correction rule to correct for errors 
only in these large variance regions. One such rule would be to 
error correct code those transform samples which correspond to 
positions in the transform domain where the transform spectrum of 
the covariance function indicates a high probability of large sample 
value* This technique alone requires an increase in bandwidth to 
facilitate the error correction. However, it has been found that the 
small increase in bandwidth in the transform domain will result in 
better reconstructions than the same increase in the spatial domain. 
Figure 7-8 demonstrates this situation where a 3.5:1 increase in 
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a . Spatial domain errors 


b. Spatial domain 
error correction 
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bandwidth has greatly improved the transform coded image over 

the spatial coded image. It is important to emphasize that the 

coding technique used for the transform domain should be tailored 

to a particular channel capacity. If the channel noise has an error 
-4 

rate less than about 10 , then it appears that no error correction 

is necessary as in Figures 7-4a, 7-5a r 7-6a, and 7-7a. However, 
under the circumstances of a high error rate, it often becomes 
desirable to transmit as many error corrected samples as 
possible at the expense of either increased bandwidth or of not trans- 
mitting the entire transform plane. Using such a system, corrected, 
but not necessarily errorless, data could be received until either all 
data (and parity bits) are received for a complete picture or until 
normal picture bandwidth has been reached, at which time trans- 
mission is terminated. In order to implement such a scheme, an 
error correcting code must be selected. 

A specific example of the potential of the transform error 

correction coding technique is presented below* A high error rate 

-2 

channel is assumed with rate Pe = 4 x 10 , Three experiments are 

implemented, one of which uses an increased bandwidth and the other 
two utilize an equal bandwidth criterion such that the exact same 
number of bits is necessary to transmit the spatial domain as the 
transform coded domain, (256) (256) (6), A Bose Chandhuri-Hocquenghem 
(BCH) code [41, p. 163] which is capable of correcting a total of 
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seven errors is 31 bits long with 6 information bits (31,6). 
Utilising an error correcting code capable of seven error correc- 
tions does not mean that the six information bits will be received 


over the noisy channel error free. Since each code word length 
has been Increased to thirty-one bits, eight or more errors per code 
word cannot be guaranteed to be corrected. The probability of 
having eight or more errors in the BCH code (31,6) is given by the 
partial sum of the binomial distribution 


where p is the binary symmetric channel error rate. This probability 
is an upper bound for the incorrect reception of a code word since the 
possibility of correct reception for greater than seven errors still 
exists but is unknown. For the specific channel error rate of 


Figure 7-8 presents the results of an experiment in which an 


increased bandwidth has been allowed to compensate for the parity 
bits necessary in the error correction code. However, a (31/61:1 
increase in bandwidth would be necessary to completely transmit the 
full data of either the space or transform domain. Allowing only a 
3, 5:1 bandwidth increase means that not all the data in the transform 


31 



p (8 or more errors) = S 
i=8 


-2 

4 X 10 the error corrected data samples will be received with 
probability of error no greater than 2, 26 x 10 ^ [42] * 
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domain is transmitted and thus Figures 7-8c and 7 - 8d are 1,44:1 

-5 

low pass sequency and frequency filters with error rate Z. 26 * 10 
The spatial domain error correction is the average of 70% BCH code 
and 30% no error coding. 

For spacecraft implementations it is desirable to transmit 

the error corrected image with no increased bandwidth requirement 

over conventional spatial domain transmission. Thus a (31/ 6) 1 1 low 

-S 

pass frequency or sequency filter with error rate 2. 26 X 10 will 
result in an equal bandwidth requirement. The results of this 
experiment are displayed in Figures 7- 9b and 7-9d, 

Because zonal low pass transform filtering is a non-adaptive 
technique for bandwidth reduction, it is desirable to utilize the adap- 
tive feature of threshold coding as a means of more optimally com- 
pensating for the parity bits necessary for error coding. Thus a run 

length coding technique utilizing 4 position bits will be used in the 
* 

transform domain. Consequently, 4 position and 6 data bits will be 
used as information bits in the transform domain for run length coded 
thresholded transform samples. Thus a new error correcting code is 
necessary and a convenient candidate is a BCH (31, 11) code. This 


A pseudo- run length coding technique is alluded to here enabling 
4 rather than 8 bits necessary for position coding [43], 
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a, Fourier run length 
error corrected re transformation 


b* Fourier zonal error 
corrected retransformation 




c* Hadamard run length 
error corrected ^transformation 


d* Hadamard zonal error 
corrected re transformation 



e* Spatial domain errors 


Figure 7-3 Surveyor Box Equal Bandwidth 

Error Correction Technique , Pe = 4 x 1 0 
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code, again of length 31 bits, has 11 information bits and is capable 
of correcting 5 or less bit errors. Consequently, the probability of 
having 6 or more errors in the BCH code for each sample is given by 
the partial sum of the binomial distribution 

3] 

p {6 or more errors) = Z) 

i~ 

“3 ? 

and is equal to 1.27 X 10 for a channel with error rate 4 X IQ - [43]. 

Thus the cost of run length coding has changed the effective error 
rate from 2* 26 x 10 to 1.27 X 10 for this example. Figure 7-9a 
and 7- 9c are the run length error corrected retransformations with 
a 5:1 bandwidth reduction to compensate for the (31/6) *1 parity 
information bandwidth increase. Consequently, again, an equal 
bandwidth criterion has been maintained. 

It is suggested that other coding techniques could be developed 
which would improve upon these results. In fact, for potential hard- 
ware systems, research ought to be undertaken to develop the best 
code for the channel error rate, bandwidth, and computational com- 
plexity allowable. 


i,. ,31-i 

P O-P) 


(7-2) 
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8. Summary 


This report has presented a theoretical development of 
several two dimensional transforms that are potentially useful for 
image coding. Of the transforms analyzed* the Fourier, Hadamard, 
and Karhunen- Loeve transforms have proven to possess the desired 
property of image energy compaction in the transform domain. 

The energy compaction property of these three transforms 
has been exploited to achieve a sample reduction by two means i 
zonal sampling and threshold sampling. In zonal sampling a sampling 
mask corresponds to the positional ordering of the largest eigenvalues 
of Idle covariance matrix of the class of images to be coded. For the 
Fourier and Hadamard transforms the best sampling mask has a 
hyperbolic shape in the transform domain. Examples of the sample 
reduction achievable by zonal sampling with the three transforms are 
shown in Figures 8- la to 8-lc. The transforms were taken in blocks 
of 16 by 16 elements. The other technique of sample deletion, called 
threshold sampling* simply entails the coding of each transform domain 
sample that exceeds a magnitude threshold level* By this technique the 
reconstruction of a particular image will suffer tiie least degradation 
from the standpoint of energy loss* Figures 8- Id to 8- If illustrate 
the performance of threshold coding. The important conclusions to 
be drawn from Figure 8-1 and the supporting experimental results 
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Fourier, hyperbolic zonal sampling 
4:1 sample reduction 


b. 



Fourier, threshold sampling 
5:1 sample reduction 



Hadamard, hyperbolic zonal sampling 
4:1 sample reduction 



♦ Karhunen-Loeve, zonal sampling 
4:1 sample reduction 



d. Hadamard, threshold sampling 
5:1 sample reduction 



, Karhunen-Loeve , threshold sampling 
5:1 sample reduction 


Figure 6-1 Summary of Fourier, Hadamard, and Karhunen-Loeve 
Transform Image Coding in 16 by 16 Element Blocks 
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of Section 3 are that; 


a* significant sample reduction factors can be obtained by zonal 
and threshold sampling in the transform domain for the 
Karhunen -Loeve, Fourier, and Hadamard transforms. 

b. threshold sampling provides better performance (higher sample 
reduction factors for the same degree of image quality) than 
zonal sampling. 

c. the Karhunen- Loeve transform exhibits somewhat better 
performance than the other two transforms, which in turn, 
exhibit about the same degree of performance. 

d. the sample reductions achieved were obtained by transform 
coding in blocks of only 16 by 16 elements. Image trans- 
formation in such small blocks can be implemented quite 
simply. 

Fast computational algorithms exist for the Fourier and 
Hadamard transforms. Computation of these transforms on a 
general purpose computer in blocks of up to 1024 by 1024 elements 
appears feasible from a computational standpoint. There is no 
fast computation algorithm for the Karhunen- Lo eve transform. This 
fact coupled with the realization that the Karhunen- Loeve does not 
perform appreciably better than the Fourier and Hadamard transforms 
seems to limit the practical utility of the Karhunen- Loeve transform. 
For these reasons the detailed analysis of the report has been limited 
primarily to the Fourier and Hadamard transforms. 
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An analysis has been performed to determine the optimum 
means of transform domain sample quantization. The results of 
this analysis indicate that for a quantization strategy in which 
each sample is coded to the same number of levels, the optimum 
quantizer places the quantization levels along a nonlinear scale 
both in sample amplitude and position in the transform domain. 
Unfortunately, the optimum is difficult to implement. Therefore, 
several nonlinear scales that could be deterministically computed 
were analyzed. The best performance has been obtained with a 
Gaussian error function quantizer. With this quantizer, good quality 
reconstructions have been obtained with 64 quantization levels 
(6 bits) per transform sample component for both the Fourier and 
Hadamard transforms. 

Zonal and threshold sampling of quantized Fourier and Hadamard 
transforms of images has been investigated in detail for a variety 
of images. The transforms have been taken in blocks of up to 256 
by 256 elements. A position coding technique for threshold sampling 
employing run length coding has been implemented and evaluated. 

Figure 8-2 illustrates the effects of threshold coded quantized Fourier 
and Hadamard transforms of images over a full frame of 256 by 256 
ele merits . 

The conclusions to be drawn from these experiments are that: 
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a* Fourier, hyperbolic zonal sampling 
4 l 1 sample reduction 


b. Fourier, threshold sampling 
5:1 sample reduction 




c. Hadamard, hyperbolic zonal sampling 
4:1 sample reduction 


d. Hadamard, threshold sampling 
5:1 sample reduction 


Figure 8-2 Summary of Fourier and Hadamard Transform Full Frame 
Image Coding — Quantized and Coded Images 
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a. for any size block, threshold sampling provides better per- 
formance than zonal sampling. 

b. performance .is better for larger size blocks, but the difference 
in performance between blocks of 16 by 16 elements and blocks 
of 256 by 256 elements is not great. 

c. for threshold sampling, simple run length coding can be 
employed to code the position of significant samples; the run 
length coding does not affect image quality, and can be 
accomplished with a relatively few number of bits per image 
element. 

The effect of channel errors on transform coded images has 
been studied. It has been found that channel errors in the transform 
domain tend to cause a small overall loss in resolution; there are 
no discrete effects like the "salt and pepper" errors that appear in 
normal spatial domain coding. Experiments verify that errors in 
the position bits coding the position of significant samples in 
threshold coding are not serious. Errors in the lowest spatial 
frequencies (sequencies) have been found to degrade an image the 
most. By applying channel error correction to a relatively small 
number of these transform domain samples, a relatively large 
improvement in the tolerance to channel errors can be obtained. 

The equivalent amount of error correction in the spatial domain 
would provide no worthwhile improvement. 
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In final summary, it can be said that Fourier and Hadamard 
transform image coding techniques are a feasible means of obtaining 
significant bandwidth compressions for digital image transmission. 
Side benefits of transmitting the Fourier or Hadamard transform of 
an image rather than the image itself are an improved tolerance to 
channel errors and the fact that image enhancement methods can be 
readily performed in the transform domain. 
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9. Recommendations' 


The general concept of transform coding has now been studied 
and evaluated rather thoroughly in this research study and by other 
investigators. There remain three areas, listed below, that merit 
further study. 

Transform Domain Coding 

Zonal sampling in the transform domain has the advantage of 
simplicity, but achievable performance is not as great as can be 
obtained by threshold sampling. However, threshold sampling requires 
position coding of significant samples. It appears that advantages of 
both techniques might be obtained by a hybrid scheme of zonal sampling* 
a set of the low spatial frequencies (sequencies) and threshold sampling 
the remainder of the transform domain. Schemes for performing this 
type of sampling should he investigated in conjunction with a study of 
the best means of position coding significant samples. 

The quantization technique presented in this report adopted the 
strategy of assigning the same number of bits per transform domain 
sample and then determining the optimum scaling of quantization levels. 
Another technique that has been reported [37] utilizes a linear quanti- 
zation scale for each sample, but the number of bits per sample is 
optimally selected to minimize the total number of image code bits for 
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a given error criterion. It appears that it would be advantageous 
to combine both strategies: assign the number of bits per sample 
on the basis of the sample variance, and select the quantization 
levels according to a nonlinear scale based upon the variance. This 
quantization method should be studied further. 

Implementation 

A number of companies have available equipment to perform 
a fast Fourier transform in one dimension for up to about 1024 
points. A few companies have built fast Hadamard transform 
devices for one dimensional transforms. There are presently no 
two dimensional transform processors on the market. 

In view of the great potential for image transform coding it 
would seem worthwhile to implement prototype Fourier and Hadamard 
transform processors. As a first step a 16 by 16 element processor 
should be built and evaluated. 

Color Image Coding 

Conventional color images are represented by three overlapping 
intensity planes corresponding to three primary colors --red, green, 
and blue. In normal television practice a linear color transition is 
made into three planes which represent the luminance (the monochro- 
matic representation of an image) and the two chrominance variations 
of the image. The spatial frequency response of the eye to chrominance 
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information is very poor. Therefore, a great deal of spatial 
low pass filtering on the chrominance 'planes can be tolerated. 
Fourier and Hadamard zonal low pass filtering appear ideal for 
this application. Studies are needed to determine the effects of 
Fourier and Hadamard filtering on the chrominance planes and to 
determine the best color transitions for the subsequent filtering 
operation. 
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APPENDIX A 


Imaere Covariance Function* 

Measurements have been made of the covariance function of an image 
to determine the fit of the Gauss -Markov process model. Figure A-l shows 
plots of the correlation between elements along a line, between elements along 
a column of the image, and between elements along the diagonal of an image. 
All measurements have been made on the Surveyor spacecraft scene. The 
data points have been fit by functions of the form A n where A is the corre- 
lation between adjacent elements and n is the separation between elements. 
The fit along the rows and columns of the image appears to be reasonably 
good. As shown in the figure there is a small deviation between the Gauss - 
Markov process model for diagonal elements and actual measurements. 


Measurements have been performed by Professor Lee D. Davisson of 
the University of Southern California. 
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CORRELATION 

ALONG A LINE 




(. 953 ) 



References 


lo W.K. Pratt U A Bibliography on Television Bandwidth Reduction 
Studies' 1 , IEEE Transactions on Information Theory , Vol. IT -13, 

No, 1 (January, 1967), pp. 114-115, 

2, A, Rosenfeld 1 'Bandwidth Reduction Bibliography”, IEEE 
Transactions on Information Theory , VoL IT- 14, No, 4 
(July, 1968), pp, 601-602. 

3. Special issue- on redundancy reduction Proceedings IEEE , Vol. 55, 
No. 3 (March, 1967). 

4. H.C, Andrews and W.K. Pratt "Fourier Transform Coding of 
Images”, Hawaii International Conference on System Sciences , 
(January, 1968), pp. 677-679* 

5, H.C, Andrews and W, K, Pratt "Television Bandwidth Reduction 

by Fourier Image Coding”, Society of Motion Picture and Television 
Engineers , 10 3rd Technical Conference , (May, 1968). 

6, H.C. Andrews and W.K. Pratt "Television Bandwidth Reduction 

by Encoding Spatial Frequencies", Journal Society of Motion Picture 
and Television Engineers , Vol. 77 (December, 1968), pp, 1279- 
1281. ~ ~~ 

7. W.K. Pratt, J. Kane, and H. C. Andrews "Hadamard Transform 
Image Coding", Proceedings IEEE , Vol. 57, No. 1 (January, 1969). 

W.K. Pratt and H.C. Andrews "Application of Fourier -Hadamard 
Transformation to Bandwidth Compression”, MIT Symposium on 
Picture Bandwidth Compression , (April, 1969). 

9. W.K. Pratt and H. C. Andrews "Two Dimensional Transform Coding 
of Images”, 1969 International Symposium of Information Theory 
Institute of Electrical and Electronic Engineers , (November, 1968). 

10. H.C. Andrews and W.K. Pratt "Transformation Coding for Noise 
Immunity and Bandwidth Reduction”, Second Annual Hawaii 
International Conference on System Sciences, (January, 1969). 

11. H.C. Andrews and W.K. Pratt "Transform Image Coding”, PIB 
International Symposium on Computer Processing in Communi - 
cations , (April, 1969). 


- 147 - 



12 . 


J.W. Cooley, P.A. W. Lewis, and P. D. Welch "Historical 
Notes on the Fast Fourier Transform", Proceedings IEEE , 

Vol. 55, October, 1967, pp. 1675-1677.’ 

13. H. C. Andrews "Fourier Coding of Images", University of 
Southern California, USCEE Report No. 271, (June, 1968). 

14. J. Hadamard "Resolution d'une Question Relative aux 
Determinants", Bulletin des Sciences Mathematiques , (2), 

Vol. 17, part 1, (1893), pp. 240-246. 

15. H.J. Ryser, Combinatorial Mathematics , John Wiley, New 
York, (1963). 

16 . S. W. Golomb, et al. Digital Communications, Prentice -Hall, 
(1964). 

17. H.F. Harmuth "A Generalized Concept of Frequency and Some 
Applications", IEEE Transactions on Information Theory , Vol. 
IT-14, No. 3, (May, 1968), pp. 375-382. 

18. J. L. Walsh "A Closed Set of Orthogonal Functions", American 
Journal Mathematics , Voi. 45, (1923), pp. 5-24. 

19. N.J, Fine "On the Walsh Functions", Transactions American 
Mathematical Society, Vol. 65, (1949), pp. 372-414. 

20. N.J. Fine "The Generalized Walsh Functions", Transactions 
American Mathematical Society , Vol. 69, (1950), pp. 66-77. 

21. G.W. Morgenthaler "On Walsh-Fourier Series", Transactions 
American Mathematical Society , Vol. 84, (T957), pp. 472-507. 

22. K.W. Henderson "Some Notes on the Walsh Functions", IEEE 
Transactions on Electronic Computers , Vol. EC-43, (February, 
1964), pp. 50-52. 

23. H. Rademacher "Einige Satze von Allgemeinen Orthogonal- _ 
Funktionen", Mathematics Annals, Vol. 87, (1922), pp. 122-138. 

24. H.C. Andrews and J. Kane "Kronecker Matrices, Computer 
Implementation, and Generalized Spectra", Journal of -the 
Association of Computer Machinery ,' (April, 1970). 


- 148 - 



25. H. C. Andrews and K. L. Caspari !r A Generalized Technique for 

Spectral Analysis", IEEE Transactions on Computers , Vol. C-9, 
No. 1 (Jaunary, 1970), pp, 16-25. ~ 

26. H.C. Andrews and W.K. Pratt "Transform Data Coding", 

PIB Symposium on Computer Proce ssing in Communications, 

April, 1969. 

27. W.K. Pratt and H. C. Andrews "Transform Processing and Coding 
of Images", University of Southern California, Electronic Sciences 
Laboratory, USCEE Report No. 341 (March, 1969), Chapter 2. 

28. A. Haar "Zur Theorie des Orthogonalen Funktionen-Systeme", 
Inaugural dissertation, Math. Annals , Vol. 69 (1910), pp. 331- 
371 and Vol. 71 (1912), pp. ‘33-53. 

29. C. Wateri "A Generalization of Haar Functions", Tohoku 
Mathematical Journal . Vol. 8 (1956), pp. 286-290. 

30. H.P. Kramer and M. V. Mathews "A Linear Coding for Trans- 
mitting a Set of Correlated Signals", IRE Transactions on 
Information Theory , Vol. IT-2 (September, 1956), pp. 41-46. 

31. J.E. Whelchel, Jr., and D. F. Guinn "The Fast Fourier-Hadamard 
Transform and Its Use in Signal Representation and Classification", 
EASCON 1968 Convention Record, (1968), pp. 561-573, 

32. J.J.Y. Huang and P.M. Schutheiss "Block Quantization of 
Correlated Gaussian Random Variables", IEEE Transactions 

on Communication Systems , Vol. CS-11, No. 3 (September, 1963), 
pp. 289-296. 

33. T.Y. Young and W. H. -Huggins "On the Representation of Electro- 
cardiographs", IEEE Transactions on Bio-Medical Electronics , 

Vol. BME-10, No. 3 (July, 1963), pp. 86-95. 

34. L.M. Goodman "A Binary Linear Transformation for Redundancy 
Reduction" , Proceedings IEEE Letters , Vol. 55, No. 3 (March 
1967), pp. 467-468. 

35. C. A. Andrews, J.M. Davies, and G. R. Schwarz "Adaptive 
Data Compression", Proceedings IEEE , Vol. 55, No. 3 
(March, 1967), pp. 267-277. 


- 149 - 



36. 


C.J. Palermo, R. V. Palermo, and H. Horowitz "The Use of 
Data Omission for Redundancy Removal", Record International 
Space Electronics and Telemetry Symposium , (1965), pp. (11)D1- 
(1 1)D16. 

37. A. Habibi and P. Wintz "Optimum linear Transformations for 
Encoding 2-Dimensional Data". 

38. D. A. Pipes Matrix Methods in Engineering, Prentice -Hall, 
Englewood Cliffs, New Jersey (1963). 

39. A. Papoulis Probability, Random Variables, and Stochastic 
Processes , McGraw-Hill Book Company (1965). 

40. P.F. Panter and W. Dite "Quantization Distortion in Pulse 
Count Modulation with Nonuniform Spacing of Levels", Pro- 
ceedings IRE, Vol. 39, No. (January, 1951), pp. 44-48. 

41. ,W. W. Peterson Error Correcting Codes ' , The MIT Press, 
Cambridge, Massachusetts (1961). 

42. "Tables of the Binomial Probability Distribution", Departments 

of Commerce, National Bureau of Standards, Applied Mathematics 
Series No. 6 (January, 1950). 

43. W.K. Pratt "Stop Scan Edge Detection Systems of Television 
B andwidth’ Reduction", USCEE Report 13IT, (June, 1965). 


-150- 




